과제의 목표는 공공 자전거의 활성화와 불만족도를 낮추는 것이다. 이 해당 목표에 대해 크게 3가지의 방향을 갖고 분석을 진행하였다.
공공성을 가지는 과제의 목표가 마음에 들었고, 공간 데이터를 기반으로 새로운 데이터 분석에 도전하고 싶었다. 도전을 통해 개인의 성장과 더불어 좋은 결과를 도출하여 공공의 이익에 일조하고 싶다.
용어 정의
총 자전거 이용량 (반입량+반출량)의 일별 평균 : 절대 수요
총 자전거 변화량 (반입량-반출량)의 일별 평균 : 상대 수요
실행 후 런타임 다시시작하세요!
!pip install -q folium==0.8.3
!pip install -q geopandas
!pip install haversine
!pip install -q mapclassify==2.3.0
# 기본 데이터 라이브러리 로드
import sys # 시스템 파라미터에 접근할 수 있게 도와준다.
print("Python version: {}". format(sys.version))
import pandas as pd # 데이터 정제에 도움을 주는 라이브러리
print("pandas version: {}". format(pd.__version__))
import matplotlib # 매트랩에서 사용하는 시각화 도구를 사용할 수 있게 도와주는 시각화 도구
print("matplotlib version: {}". format(matplotlib.__version__))
import numpy as np # 행렬 계산을 위해 필요한 라이브러리
print("NumPy version: {}". format(np.__version__))
import scipy as sp # 수학 관련 함수가 내장된 라이브러리
print("SciPy version: {}". format(sp.__version__))
import IPython
from IPython import display # 주피터 노트북에서 예쁘게 시각화 해주는 도구
print("IPython version: {}". format(IPython.__version__))
import sklearn # 각종 통계 도구와 머신 러닝 알고리즘이 내장되어 있는 라이브러리
print("scikit-learn version: {}". format(sklearn.__version__))
import tensorflow as tf
print("tensorflow version: {}".format(tf.__version__))
import geopandas as gpd # 공간정보 라이브러리
print("geopandas version: {}".format(gpd.__version__))
# 파이썬 내장 라이브러리
import random
import datetime
import time
import os
from datetime import datetime, timedelta
from dateutil.relativedelta import relativedelta
import copy
# 모델링 라이브러리, 수학 계산 툴 로드
# 일반적인 모델링 라이브러리
from sklearn import svm, tree, linear_model, neighbors, naive_bayes, ensemble, discriminant_analysis, gaussian_process
from xgboost import XGBClassifier
# 모델링 시 헬퍼 함수들
from sklearn.preprocessing import OneHotEncoder, LabelEncoder
from sklearn import feature_selection
from sklearn import model_selection
from sklearn import metrics
# interpolation 함수
import scipy.interpolate as spi
from scipy import stats
# model load
import pickle
# from sklearn.externals import joblib
# 시각화 도구
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
import seaborn as sns
from pandas.plotting import scatter_matrix
# 시각화 도구 default 세팅
# 주피터 노트북에서 plot 결과를 볼 수 있게 해준다.
%matplotlib inline
mpl.style.use('ggplot') # matplotlib에서 plot되는 결과를 선택할 수 있다.
sns.set_style('white') # seaborn에서 사용할 style을 설정할 수 있다.
pylab.rcParams['figure.figsize'] = 12,8 # plot의 크기와 선 등의 기본 값을 설정할 수 있다.
plt.style.use(['fivethirtyeight'])
sns.set_style('darkgrid')
from IPython.display import display #print가 아닌 display()로 연속 출력
from IPython.display import HTML #출력 결과를 HTML로 생성
import fiona #공간데이터를 딕셔너리 형태 등으로 접근할 수 있는 라이브러리
# 좌표계 정의/변환용 라이브러리
import pyproj
# 좌표간의 거리 구하기 위한 라이브러리
from haversine import haversine
# Jupyter Notebook 이나 ipython 을 사용하다보면 향후 버전이 올라갈 때 변경될 사항 등을 알려주는 경고 메시지(warning message)를 뜨지 않게 해준다.
import warnings
warnings.filterwarnings('ignore')
print('-'*25)
# 뒤에서 overlay를 활용하기 위함. 설치 후 런타임 다시 시작할 것
import rtree
import pathlib
from geoband import API
import mapclassify
import os
# 사용하는 디렉토리 정의
BASE_WORKING_DIR = os.getcwd()
DATA_RAW_PATH = os.path.join(BASE_WORKING_DIR, "input")
# 사용 좌표계
default_crs = "epsg:4326"
meter_crs = "epsg:5179"
# 데이터 불러오기
input_path = pathlib.Path('./input')
if not input_path.is_dir():
input_path.mkdir()
API.GetCompasData('SBJ_2007_001', '1', input_path.joinpath('운영이력.csv'))
API.GetCompasData('SBJ_2007_001', '2', input_path.joinpath('자전거스테이션.csv'))
API.GetCompasData('SBJ_2007_001', '3', input_path.joinpath('꽃박람회일정.csv'))
API.GetCompasData('SBJ_2007_001', '4', input_path.joinpath('KINTEX행사일정.csv'))
API.GetCompasData('SBJ_2007_001', '5', input_path.joinpath('기상정보.csv'))
API.GetCompasData('SBJ_2007_001', '6', input_path.joinpath('인구(거주)분포도(100M X 100M).geojson'))
API.GetCompasData('SBJ_2007_001', '7', input_path.joinpath('인구통계.csv'))
API.GetCompasData('SBJ_2007_001', '8', input_path.joinpath('행정경계(시군구).geojson'))
API.GetCompasData('SBJ_2007_001', '9', input_path.joinpath('행정경계(읍면동).geojson'))
API.GetCompasData('SBJ_2007_001', '10', input_path.joinpath('도시계획(공간시설).geojson'))
API.GetCompasData('SBJ_2007_001', '11', input_path.joinpath('도시계획(공공문화제육시설).geojson'))
API.GetCompasData('SBJ_2007_001', '12', input_path.joinpath('도시계획(교통시설).geojson'))
API.GetCompasData('SBJ_2007_001', '13', input_path.joinpath('용도지역지구(습지보호지역).geojson'))
API.GetCompasData('SBJ_2007_001', '14', input_path.joinpath('고양시 지적도.geojson'))
API.GetCompasData('SBJ_2007_001', '15', input_path.joinpath('도로명주소_건물.geojson'))
API.GetCompasData('SBJ_2007_001', '16', input_path.joinpath('도로명주소_도로.geojson'))
API.GetCompasData('SBJ_2007_001', '17', input_path.joinpath('일반건물 분포도(100M X 100M).geojson'))
API.GetCompasData('SBJ_2007_001', '18', input_path.joinpath('행사장_공간정보.csv'))
API.GetCompasData('SBJ_2007_001', '19', input_path.joinpath('전철역_공간정보.csv'))
API.GetCompasData('SBJ_2007_001', '20', input_path.joinpath('고양시 버스정류소.csv'))
API.GetCompasData('SBJ_2007_001', '21', input_path.joinpath('버스 정류장별 승하차 정보.csv'))
API.GetCompasData('SBJ_2007_001', '22', input_path.joinpath('주차장정보.csv'))
API.GetCompasData('SBJ_2007_001', '23', input_path.joinpath('고양시덕양구_DEM.img'))
API.GetCompasData('SBJ_2007_001', '24', input_path.joinpath('고양시일산동구_DEM.img'))
API.GetCompasData('SBJ_2007_001', '25', input_path.joinpath('고양시일산서구_DEM.img'))
API.GetCompasData('SBJ_2007_001', '26', input_path.joinpath('고양시 공연장 박물관 정보.csv'))
API.GetCompasData('SBJ_2007_001', '27', input_path.joinpath('고양시 체육시설 현황 정보.csv'))
API.GetCompasData('SBJ_2007_001', '28', input_path.joinpath('코드정의서.xlsx'))
API.GetCompasData('SBJ_2007_001', '29', input_path.joinpath('지하철 역별 이용객수.csv'))
API.GetCompasData('SBJ_2007_001', '30', input_path.joinpath('고양시_덕양구_고도.geojson'))
API.GetCompasData('SBJ_2007_001', '31', input_path.joinpath('고양시_일산동구_고도.geojson'))
API.GetCompasData('SBJ_2007_001', '32', input_path.joinpath('고양시_일산서구_고도.geojson'))
API.GetCompasData('SBJ_2007_001', '33', input_path.joinpath('고양시_인도.geojson'))
API.GetCompasData('SBJ_2007_001', '34', input_path.joinpath('행정경계(행정동기준).geojson'))
API.GetCompasData('SBJ_2007_001', '35', input_path.joinpath('고양시_도시화지역경계.geojson'))
for path in list(input_path.glob('*.csv')) + list(input_path.glob('*.geojson')):
print(path)
각 지역과 도로에 대한 코드의 의미를 정리한 테이블이다.
| code definition | dataframe variable name |
|---|---|
| 읍면동 코드(EMD_CD) | dong_code_df |
| 시군구 코드(SIG_CD) | gu_code_df |
| 광역도로 구분 코드(WDR_RD_CD) | road_code1_df |
| 도로구간종속구분코드(RDS_DPN_SE) | road_code2_df |
| 도로위계구분코드(ROA_CLS_SE) | road_code3_df |
| 건축물용도코드(BDTYP_CD) | building_usage_code_df |
| code definition | dataframe variable name |
|---|---|
| 법정동 코드 -> 이름 | gu_legal_dic |
| 행정동 코드 -> 이름 | gu_administration_dic |
| 법정동 시군구 code | value |
|---|---|
| '41281' | 덕양구 |
| '41285' | 일산동구 |
| '41287' | 일산서구 |
| 행정동 시군구 code | value |
|---|---|
| '31100' | 고양시 |
| '31101' | 덕양구 |
| '31103' | 일산동구 |
| '31104' | 일산서구 |
# 코드 정의서 로드
code_def_df = pd.read_excel(os.path.join(DATA_RAW_PATH, "코드정의서.xlsx"))
# 읍면동 코드(EMD_CD)
dong_code_df = code_def_df.loc[2:55]
dong_code_df.columns = ["code", "meaning"]
dong_code_df = dong_code_df.reset_index(drop=True)
dong_code_df.head()
# 시군구 코드(SIG_CD)
gu_code_df = code_def_df.loc[553:555]
gu_code_df.columns = ["code", "meaning"]
gu_code_df = gu_code_df.reset_index(drop=True)
gu_code_df.head()
# 구코드 법정동 기반 dictionary를 제작해둔다.
gu_code_list = list(gu_code_df["code"].apply(lambda x: str(x)))
gu_mean_list = list(gu_code_df["meaning"].apply(lambda x: x.split()[2]))
gu_legal_dic = dict(zip(gu_code_list, gu_mean_list))
gu_legal_dic
# 구코드 행정동 기반 dictionary를 제작해둔다.
gu_code_list = ["31100", "31101", "31103", "31104"]
gu_mean_list = ["고양시", "덕양구", "일산동구", "일산서구"]
gu_administration_dic = dict(zip(gu_code_list, gu_mean_list))
gu_administration_dic
# 광역도로 구분 코드(WDR_RD_CD)
road_code1_df = code_def_df.loc[58:60]
road_code1_df.columns = ["code", "meaning"]
road_code1_df = road_code1_df.reset_index(drop=True)
road_code1_df.head()
# 도로구간종속구분코드(RDS_DPN_SE)
road_code2_df = code_def_df.loc[540:542]
road_code2_df.columns = ["code", "meaning"]
road_code2_df = road_code2_df.reset_index(drop=True)
road_code2_df.head()
# 도로위계구분코드(ROA_CLS_SE)
road_code3_df = code_def_df.loc[546:549]
road_code3_df.columns = ["code", "meaning"]
road_code3_df = road_code3_df.reset_index(drop=True)
road_code3_df.head()
# 건축물용도코드(BDTYP_CD)
building_usage_code_df = code_def_df.loc[64:536]
building_usage_code_df.columns = ["code", "meaning"]
building_usage_code_df = building_usage_code_df.reset_index(drop=True)
building_usage_code_df.head()
주어진 데이터를 기반으로 아래와 같은 geopandas dataframe 형태의 변수를 가지고 있는다.
| code definition | dataframe variable name |
|---|---|
| 정류장 중심 정리 데이터 | station_df |
| STATION_ID | NAME | LAT | LON | geometry |
|---|---|---|---|---|
| 101 | xx아파트 | 127 | 36 | POINT(x, y) |
# 데이터 불러오기
station_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "자전거스테이션.csv"))
# column명 변경
station_df.columns = ['STATION_ID', 'STATION_NM', 'HOLDER_CAPACITY', 'LAT', 'LON']
station_df.head()
# LAT, LON 를 geometry point 객체로 변환하고 변수로 추가한다.
import pyproj
from fiona.crs import from_epsg
coulumns = ['STATION_ID', 'STATION_NM', 'HOLDER_CAPACITY', 'LAT', 'LON']
station_df = gpd.GeoDataFrame(station_df[coulumns], geometry=gpd.points_from_xy(station_df.LON, station_df.LAT))
station_df.crs = default_crs
station_df.head()
주어진 데이터를 기반으로 아래와 같은 geopandas dataframe 형태의 변수를 가지고 있는다.
| 변수명 | 의미 |
|---|---|
| LEAS_NO | 대여 번호 |
| LEAS_DATE | 대여 시간 |
| LEAS_STATION | 대여 스테이션 번호 |
| LEAS_DEF_NO | 대여 거치대 번호 |
| RTN_DATE | 반납 시간 |
| RTN_STATION | 반납 스테이션 번호 |
| RTN_DEF_NO | 반납 거치대 번호 |
| TRNV_QTY | 추정 이동 거리(m) |
| MEMB_DIV | 회원 구분(비회원은 99이며 나머지는 정회원) |
| MEMB_NO | 회원 번호 |
| TEMP_MEMB_NO | 비회원 번호 |
| BIKE_TAG | 자전거 번호 |
| RTN_PROCESS | 관제반납구분값 (01:관제반납/02:관리자반납/03:관리자반출/04:타사용자 반납/05:다른 사용자 반출) |
| code definition | dataframe variable name |
|---|---|
| 운행 중심 정리 데이터 | route_df |
| LEAS_NO | LEAS_DATE | LEAS_STATION | LEAS_DEF_NO | LEAS_LAT | LEAS_LON | LEAS_geometry | LEAS_HOLDER_CAPACITY | LEAS_STATION_NAME | RTN_DATE | RTN_STATION | RTN_DEF_NO | RTN_LAT | RTN_LON | RTN_geometry | RTN_HOLDER_CAPACITY | RTN_STATION_NAME | TRNV_QTY | MEMB_DIV | MEMB_NO | TEMP_MEMB_NO | BIKE_TAG | RTN_PROCESS | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 15945541 | 2017-01-01 00:00:41 | 213 | 18 | 37.659119 | 126.787998 | POINT (126.78800 37.65912) | 25.0 | 마두1동 사거리 | 2017-01-01 | 00:13:52 | 260 | 17 | 37.662030 | 126.768952 | POINT (126.76895 37.66203) | 30.0 | ★장항 제1공영주차장 | 0.0 | 6 | 164203 | 0.0 | 1A844000000BB7 | NaN |
# 데이터 불러오기
route_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "운영이력.csv"))
print(route_df.shape)
route_df.head()
# 이상한 변수 LEAS_STAT이 있다. 삭제해주자.
route_df = route_df.drop(columns=["LEAS_STAT"])
# 회원 비회원은 int와 string이 섞여있다.
route_df["MEMB_DIV"].unique()
# str로 형변환 해준다.
route_df["MEMB_DIV"] = route_df["MEMB_DIV"].apply(lambda x: str(x))
# LEAS_STATION, RTN_STATION에 대한 정류장 위치를 추가해주자.
change_list = ["LAT", "LON", "geometry", "HOLDER_CAPACITY", "STATION_NM"]
route_df = pd.merge(route_df, station_df, how="left", left_on="LEAS_STATION", right_on="STATION_ID")
route_df.rename(columns=lambda x: f"LEAS_{x}" if x in change_list else x, inplace=True)
route_df = route_df.drop(columns=["STATION_ID"])
route_df = pd.merge(route_df, station_df, how="left", left_on="RTN_STATION", right_on="STATION_ID")
route_df.rename(columns=lambda x: f"RTN_{x}" if x in change_list else x, inplace=True)
route_df = route_df.drop(columns=["STATION_ID"])
route_df = route_df[
[
"LEAS_NO",
"LEAS_DATE",
"LEAS_STATION",
"LEAS_DEF_NO",
"LEAS_LAT",
"LEAS_LON",
"LEAS_geometry",
"LEAS_HOLDER_CAPACITY",
"LEAS_STATION_NM",
"RTN_DATE",
"RTN_STATION",
"RTN_DEF_NO",
"RTN_LAT",
"RTN_LON",
"RTN_geometry",
"RTN_HOLDER_CAPACITY",
"RTN_STATION_NM",
"TRNV_QTY",
"MEMB_DIV",
"MEMB_NO",
"TEMP_MEMB_NO",
"BIKE_TAG",
"RTN_PROCESS",
]
]
route_df.head()
지적도에 있는 PNU를 쪼개고, 이를 기반으로 feature를 추가한 데이터 구조를 제작한다.


| code definition | dataframe variable name |
|---|---|
| 지적도 | goyang_geo_gdf |
| 구별 geometry | administration_boundary_city_gdf |
| 법정동 : 구별/동별 geometry | legal_boundary_town_gdf |
| 행정동 : 구별/동별 geometry | administration_boundary_town_gdf |
# 지역 데이터 로드
goyang_geo_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '고양시 지적도.geojson'), encoding='cp949')
administration_boundary_city_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '행정경계(시군구).geojson'), encoding='cp949')
legal_boundary_town_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '행정경계(읍면동).geojson'), encoding='cp949')
administration_boundary_town_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '행정경계(행정동기준).geojson'), encoding='cp949')
# 좌표계 설정
goyang_geo_gdf = goyang_geo_gdf.set_crs(default_crs)
administration_boundary_city_gdf = administration_boundary_city_gdf.set_crs(default_crs)
legal_boundary_town_gdf = legal_boundary_town_gdf.set_crs(default_crs)
administration_boundary_town_gdf = administration_boundary_town_gdf.set_crs(default_crs)
# 좌표계 확인
print("goyang_geo_gdf 현재 좌표계 :", goyang_geo_gdf.crs)
print("administration_boundary_city_gdf 현재 좌표계 :", administration_boundary_city_gdf.crs)
print("legal_boundary_town_gdf 현재 좌표계 :", legal_boundary_town_gdf.crs)
print("administration_boundary_town_gdf 현재 좌표계 :", administration_boundary_town_gdf.crs)
# 지적도 데이터 살펴보기
print(goyang_geo_gdf.shape)
goyang_geo_gdf.head()
# plot 해보니 outlier가 검출된다. 위도가 지나치게 높은 녀석들이 존재했다. 이를 제거하자.
# 이미 알아냈기 때문에 코드 수행할 필요가 없다.
# x = list(range(goyang_geo_gdf.shape[0]))
# y = []
# for i in range(goyang_geo_gdf.shape[0]):
# temp = 0
# for j in range(4):
# temp += np.float(goyang_geo_gdf.loc[i, "geometry"].boundary[0].coords[j][1])
# if temp/5 > 30.4:
# print(i)
# y.append(temp/5)
# # plt.plot(x, y, "o")
# 검출된 아웃라이어를 제거한다.
outlierList = [
15096,
38110,
45014,
45015,
68191,
114263,
114266,
150993,
150994,
150995,
155565,
155566,
176546,
]
goyang_geo_gdf = goyang_geo_gdf.drop(index=outlierList)
goyang_geo_gdf.plot()
# 시군구 행정경계 데이터
administration_boundary_city_gdf.head()
# geometry 완결성 검사
print("유효하지 않은 데이터 갯수 : ", administration_boundary_city_gdf.shape[0] - administration_boundary_city_gdf.geometry.is_valid.sum())
# 데이터를 plot해본다.
ax = administration_boundary_city_gdf.plot(column="SIG_CD", cmap='viridis')
ax.set_title("administration city boundary", fontsize=20)
ax.set_axis_off()
plt.show()
# 법정동 행정경계
legal_boundary_town_gdf.head()
# Code를 기반으로 구 정보를 추가한다.
legal_boundary_town_gdf["GU_CD"] = legal_boundary_town_gdf["EMD_CD"].apply(lambda x: str(x)[:5])
legal_boundary_town_gdf["GU_NM"] = legal_boundary_town_gdf["GU_CD"].apply(lambda x: gu_legal_dic[f"{x}"])
legal_boundary_town_gdf.columns = ["DONG_CD", "DONG_NM", "geometry", "GU_CD", "GU_NM"]
legal_boundary_town_gdf = legal_boundary_town_gdf[["GU_CD", "GU_NM", "DONG_CD", "DONG_NM", "geometry"]]
legal_boundary_town_gdf.head()
# geometry 완결성 검사
print("유효하지 않은 데이터 갯수 : ", legal_boundary_town_gdf.shape[0] - legal_boundary_town_gdf.geometry.is_valid.sum())
# 데이터를 plot해본다.
ax = legal_boundary_town_gdf.plot(column="DONG_NM", cmap='plasma')
ax.set_title("legal town boundary", fontsize=20)
ax.set_axis_off()
plt.show()
# 추후에 dissolve를 사용하여 통계량을 검증할 목적으로 테스트한다.
si_gdf = legal_boundary_town_gdf.dissolve(by='GU_NM').reset_index(drop=True)
si_gdf.head()
# Key는 CD로, 구역을 나누는 것은 NM으로 진행한다.
ax = si_gdf.plot(column="GU_CD", cmap='viridis')
ax.set_title("administration city boundary (legal) by dissolve", fontsize=20)
ax.set_axis_off()
plt.show()
# 행정동 행정경계
administration_boundary_town_gdf.head()
# Code를 기반으로 구 정보를 추가한다.
administration_boundary_town_gdf["GU_CD"] = administration_boundary_town_gdf["행정동코드"].apply(lambda x: str(x)[:5])
administration_boundary_town_gdf["GU_NM"] = administration_boundary_town_gdf["GU_CD"].apply(lambda x: gu_administration_dic[f"{x}"])
administration_boundary_town_gdf.columns = ["DONG_CD", "DONG_NM", "geometry", "GU_CD", "GU_NM"]
administration_boundary_town_gdf = administration_boundary_town_gdf[["GU_CD", "GU_NM", "DONG_CD", "DONG_NM", "geometry"]]
administration_boundary_town_gdf.head()
# geometry 완결성 검사
print("잘못된 데이터 갯수 : ", administration_boundary_town_gdf.shape[0] - administration_boundary_town_gdf.geometry.is_valid.sum())
# 오류가 있는 데이터를 검색한다.
administration_boundary_town_gdf.geometry[~administration_boundary_town_gdf.geometry.is_valid]
# 데이터를 plot해본다.
ax = administration_boundary_town_gdf.plot(column="DONG_NM", cmap='plasma')
ax.set_title("administration boundary", fontsize=20)
ax.set_axis_off()
plt.show()
살짝 데이터 토폴로지에 문제가 있는 것으로 보이나, 직접적으로 해당 데이터를 사용할 것이 아니기에 추후 처리는 하지 않았다.
인구 거주 분포도를 로드하고, 각 구역의 인구 수를 알기 쉽게 시각화 한다.
| code definition | dataframe variable name |
|---|---|
| 인구(거주)분포도(100M_X_100M) | population_dist_gdf |
| 월간 인구통계 | population_statistics_df |
| 동기준 인구통계 | population_statistics_dong_df |
# 데이터 로드
population_dist_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '인구(거주)분포도(100M X 100M).geojson'), encoding='cp949')
population_statistics_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "인구통계.csv"))
# 좌표계 변경
population_dist_gdf = population_dist_gdf.set_crs(default_crs)
print("population_dist_gdf 현재 좌표계 :", population_dist_gdf.crs)
population_dist_gdf.head()
# Nan은 거주인구가 0인 곳이다.
population_dist_gdf["val"] = population_dist_gdf["val"].fillna(0)
# 거주인구가 없는 곳과 있는 곳으로 나누어 시각화를 진행한다.
population_dist_not_zero_gdf = population_dist_gdf[population_dist_gdf["val"] > 0]
population_dist_zero_gdf = population_dist_gdf[population_dist_gdf["val"] == 0]
ax = population_dist_not_zero_gdf.plot( column="val", cmap='copper', scheme='Fisher_Jenks', legend=True, legend_kwds={'loc': 'lower left'}, k=10)
population_dist_zero_gdf.plot(color="silver", ax=ax)
station_df.plot(ax=ax, marker=",", color="k", markersize=10)
ax.set_title("Living Population distribution", fontsize=20)
ax.set_axis_off()
plt.show()
# 인구 통계 datatype 전처리
population_statistics_df["총인구수"] = population_statistics_df["총인구수"].apply(lambda x: int(x.replace(',' , '')))
population_statistics_df["세대수"] = population_statistics_df["세대수"].apply(lambda x: int(x.replace(',' , '')))
population_statistics_df["남자 인구수"] = population_statistics_df["남자 인구수"].apply(lambda x: int(x.replace(',' , '')))
population_statistics_df["여자 인구수"] = population_statistics_df["여자 인구수"].apply(lambda x: int(x.replace(',' , '')))
# 변수명 변경
population_statistics_df.rename(columns = {"행정구역": "EMD_KOR_NM", "조회기준": "CHECK_DATE", "총인구수": "TOTAL_POP", "세대수": "TOTAL_HOUSEHOLDS", "세대당 인구": "POP_PER_HOUSEHOLDS", "남자 인구수": "M_POP", "여자 인구수": "W_POP", "남여 비율": "MW_RATE"}, inplace = True)
population_statistics_df.head()
빌딩 연면적 분포도를 로드하고, 각 구역의 연면적을 알기 쉽게 시각화 한다.
| code definition | dataframe variable name |
|---|---|
| 일반건물 분포도(100M X 100M) | building_area_gdf |
# 건물 분포도 로드
building_area_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '일반건물 분포도(100M X 100M).geojson'), encoding='cp949')
# 좌표계 변경
building_area_gdf = building_area_gdf.set_crs(default_crs)
print("building_area_gdf 현재 좌표계 :", building_area_gdf.crs)
building_area_gdf.head()
# Nan은 연면적이 0인 곳이다.
building_area_gdf["val"] = building_area_gdf["val"].fillna(0)
building_area_gdf.head()
# 빌딩 면적 분포
building_area_not_zero_gdf = building_area_gdf[building_area_gdf["val"] > 0]
building_area_zero_gdf = building_area_gdf[building_area_gdf["val"] == 0]
ax = building_area_not_zero_gdf.plot( column="val", cmap='copper', scheme='Fisher_Jenks', legend=True, legend_kwds={'loc': 'lower left'}, k=10)
building_area_zero_gdf.plot(color="silver", ax=ax)
station_df.plot(ax=ax, marker=",", color="k", markersize=10)
ax.set_title("Building Area Distribution", fontsize=20)
ax.set_axis_off()
plt.show()
주어진 데이터를 기반으로 아래와 같은 table을 만든다.
| code definition | dataframe variable name |
|---|---|
| 행사장 공간정보 | eventhall_space_df |
| 공연장 공간정보 | concerthall_space_df |
| 박물관/미술관 공간정보 | museum_space_df |
| 버스정류장 정보 | busstop_df |
| 지하철 역 정보 | subway_space_df |
| 공공시설 정보 | city_plan_public_gdf |
| 체육시설 정보 | sport_complex_info_df |
| 주차장 정보 | parking_df |
| NAME | LAT | LON | geometry |
|---|---|---|---|
| KINTEX | 37.669028 | 126.746184 | POINT (126.74618 37.66903) |
# 행사장 로드
eventhall_space_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "행사장_공간정보.csv"))
culture_complex_info_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "고양시 공연장 박물관 정보.csv"))
# column명 변경
eventhall_space_df.columns = ['NAME', 'LAT', 'LON']
eventhall_space_df.head()
# LAT, LON 를 geometry point 객체로 변환하고 변수로 추가한다.
import pyproj
from fiona.crs import from_epsg
coulumns = ['NAME', 'LAT', 'LON']
eventhall_space_df = gpd.GeoDataFrame(eventhall_space_df[coulumns], geometry=gpd.points_from_xy(eventhall_space_df.LON, eventhall_space_df.LAT))
eventhall_space_df.crs = default_crs
eventhall_space_df.head()
culture_complex_info_df.sample(5)
# 공연장 정보만 새로운 테이블로 관리한다.
concerthall_space_df = culture_complex_info_df[culture_complex_info_df["servicetype"] == "공연장"]
concerthall_space_df.columns = ["servicetype", "NAME", "LON", "LAT"]
concerthall_space_df.head()
# LAT, LON 를 geometry point 객체로 변환하고 변수로 추가한다.
import pyproj
from fiona.crs import from_epsg
coulumns = ['NAME', 'LAT', 'LON']
concerthall_space_df = gpd.GeoDataFrame(concerthall_space_df[coulumns], geometry=gpd.points_from_xy(concerthall_space_df.LON, concerthall_space_df.LAT))
concerthall_space_df.crs = default_crs
concerthall_space_df.head()
# 박물관/미술관 정보만 새로운 테이블로 관리한다.
museum_space_df = culture_complex_info_df[culture_complex_info_df["servicetype"] == "박물관/미술관"]
museum_space_df.columns = ["servicetype", "NAME", "LON", "LAT"]
museum_space_df.head()
# LAT, LON 를 geometry point 객체로 변환하고 변수로 추가한다.
import pyproj
from fiona.crs import from_epsg
coulumns = ['NAME', 'LAT', 'LON']
museum_space_df = gpd.GeoDataFrame(museum_space_df[coulumns], geometry=gpd.points_from_xy(museum_space_df.LON, museum_space_df.LAT))
museum_space_df.crs = default_crs
museum_space_df.head()
| STATION_ID | STATION_NM | LAT | LON | geometry | GETON_CNT |
|---|---|---|---|---|---|
| 219000619 | 신일비즈니스고등학교 | 37.683296 | 126.760592 | POINT (126.76059 37.68330) | 3496.0 |
# 버스 정류장 로드
busstop_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "고양시 버스정류소.csv"))
# 최근 1년치 고양시 버스 정류장별 승하차 정보이다.
busstop_getin_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "버스 정류장별 승하차 정보.csv"))
busstop_df.columns = ["STATION_NM", "STATION_ID", "LON", "LAT"]
busstop_df = busstop_df[["STATION_ID", "STATION_NM", "LAT", "LON"]]
print(busstop_df.shape)
busstop_df.head()
# LAT, LON 를 geometry point 객체로 변환하고 변수로 추가한다.
import pyproj
from fiona.crs import from_epsg
coulumns = ["STATION_ID", "STATION_NM", "LAT", "LON"]
busstop_df = gpd.GeoDataFrame(busstop_df[coulumns], geometry=gpd.points_from_xy(busstop_df.LON, busstop_df.LAT))
busstop_df.crs = default_crs
busstop_df.head()
# 승하차 데이터 보기
print(busstop_getin_df.shape)
busstop_getin_df.head()
# 정류장 데이터에 합치기
busstop_getin_df = busstop_getin_df.groupby(['STATION_ID', 'STATION_NM'])['GETON_CNT'].sum()
busstop_getin_df = busstop_getin_df.reset_index()
print(busstop_getin_df.shape)
busstop_getin_df.head()
# 버스 정류장 데이터에 승하차 데이터를 추가한다.
busstop_df = pd.merge(busstop_df, busstop_getin_df, on = ["STATION_ID", "STATION_NM"], how = "left")
print(busstop_df.shape)
busstop_df.head()
| RAIL_NM | STATION_NM | LOT_NUM_ADDR | ROAD_NM_ADDR | LON | LAT | geometry | IN | OUT | TOTAL | NET_TOTAL |
|---|---|---|---|---|---|---|---|---|---|---|
| 3호선 | 대화 | 경기도 고양시 일산서구 대화동 2221 | 경기도 고양시 일산서구 중앙로 지하 1569 | 37.675846 | 126.747206 | POINT (37.676 126.747) | 5547363.0 | 4460662.0 | 10008025.0 | 1086701.0 |
# 지하철 역 로드
subway_space_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "전철역_공간정보.csv"))
subway_getin_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "지하철 역별 이용객수.csv"))
# 변수명 변경
subway_space_df.columns = ["RAIL_NM", "STATION_NM", "LOT_NUM_ADDR", "ROAD_NM_ADDR", "LON", "LAT"]
subway_space_df.head()
# LAT, LON 를 geometry point 객체로 변환하고 변수로 추가한다.
import pyproj
from fiona.crs import from_epsg
coulumns = ["RAIL_NM", "STATION_NM", "LOT_NUM_ADDR", "ROAD_NM_ADDR", "LAT", "LON"]
subway_space_df = gpd.GeoDataFrame(subway_space_df[coulumns], geometry=gpd.points_from_xy(subway_space_df.LON, subway_space_df.LAT))
subway_space_df = subway_space_df.set_crs(default_crs)
print(subway_space_df.shape)
subway_space_df.head()
# 대곡역 2개 중 3호선 라인을 지웠다.
subway_space_df = subway_space_df.drop([0])
subway_space_df = subway_space_df.reset_index(drop=True)
print(subway_space_df.shape)
subway_space_df.head()
# 화전역은 한국 항공대로 표현되어 있다.
subway_space_df[subway_space_df["STATION_NM"]== "화전(한국항공대)"]
# 화전(한국항공대)을 화전으로 바꾼다.
subway_space_df.loc[9, "STATION_NM"] = "화전"
subway_space_df[subway_space_df["STATION_NM"]== "화전"]
# 승하차 데이터의 전체 합계만 가져온다.
subway_getin_df = subway_getin_df.iloc[:, [1,4,12]]
subway_getin_df.head()
# 지축역은 분리해서 관리했나보다. 이 부분을 병합할 필요성이 있다.
subway_getin_df[subway_getin_df["역명"] == "지축"]
# 지하철 역에 관해 중복 제거한 승하차 데이터를 만들기 위한 빈 테이블을 만든다.
subway_getin_temp_df = pd.DataFrame(np.zeros((int(len(subway_getin_df)/3), 5)), columns= ["STATION_NM", "IN", "OUT", "TOTAL", "NET_TOTAL"])
len(subway_getin_temp_df)
# 지하철 역에 관해 중복을 제거한 승하차 데이터를 만든다.
for i in range(len(subway_getin_temp_df)):
subway_getin_temp_df["STATION_NM"][i] = subway_getin_df.iloc[i*3,0]
subway_getin_temp_df["IN"][i] = subway_getin_df.iloc[i*3,2]
subway_getin_temp_df["OUT"][i] = subway_getin_df.iloc[i*3+1,2]
subway_getin_temp_df["TOTAL"][i] = subway_getin_df.iloc[i*3 + 2,2]
subway_getin_temp_df["NET_TOTAL"][i] = (subway_getin_df.iloc[i*3,2] - subway_getin_df.iloc[i*3 + 1,2])
subway_getin_temp_df.loc[21] = subway_getin_temp_df[subway_getin_temp_df["STATION_NM"] == "지축"].sum()
subway_getin_temp_df.iloc[21,0] = "지축역"
subway_getin_temp_df = subway_getin_temp_df.drop(index=[0,11])
subway_getin_temp_df = subway_getin_temp_df.reset_index(drop=True)
print(subway_getin_temp_df.shape)
subway_getin_temp_df.head()
# 최종적으로 승하차 데이터를 지하철 위치 데이터에 합한다.
subway_space_df = pd.merge(subway_space_df, subway_getin_temp_df, how="left")
subway_space_df.head()
| TYPE | AREA | LENGTH | geometry | |
|---|---|---|---|---|
| 기타 | 사회복지시설 | 2627 | 211 | MULTIPOLYGON (((126.86736 37.62940, 126.86761 ... |
# 공공시설 로드
city_plan_public_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '도시계획(공공문화제육시설).geojson'), encoding='cp949')
# 좌표계 변경
city_plan_public_gdf = city_plan_public_gdf.set_crs(default_crs)
print("city_plan_public_gdf 현재 좌표계 :", city_plan_public_gdf.crs)
# 데이터 확인
print(city_plan_public_gdf.shape)
city_plan_public_gdf.head()
# 변수명 변경
city_plan_public_gdf.columns = ["TYPE", "AREA", "LENGTH", 'geometry']
city_plan_public_gdf.head()
# 공공 시설에 해당하는 유형은 다음과 같다
city_plan_public_gdf["TYPE"].unique()
# 이를 간단하게 plot한다. 골고루 분포되어 있다.
city_plan_public_gdf.plot()
| TYPE | NAME | LAT | LON | geometry |
|---|---|---|---|---|
| 육상 | 고양종합운동장(주,보조) | 37.675787 | 126.742594 | POINT (126.74259 37.67579) |
# 체육시설을 로드한다.
sport_complex_info_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "고양시 체육시설 현황 정보.csv"))
# 필요한 변수만을 가져온다.
sport_complex_info_df = sport_complex_info_df[["GBN", "FAC_NAME", "X", "Y"]]
sport_complex_info_df.head()
# 변수명을 변경한다.
sport_complex_info_df.columns = ["TYPE", "NAME", "LON", "LAT"]
sport_complex_info_df.head()
# LAT, LON 를 geometry point 객체로 변환하고 변수로 추가한다.
import pyproj
from fiona.crs import from_epsg
coulumns = ["TYPE", "NAME", "LAT", "LON"]
sport_complex_info_df = gpd.GeoDataFrame(sport_complex_info_df[coulumns], geometry=gpd.points_from_xy(sport_complex_info_df.LON, sport_complex_info_df.LAT))
sport_complex_info_df.crs = default_crs
sport_complex_info_df.head()
| NAME | ROAD_NM_ADDR | LOT_NUM_ADDR | SPACE | LON | LAT | geometry |
|---|---|---|---|---|---|---|
| 고양시청 부설주차장 | 경기도 고양시 덕양구 고양시청로10 | 경기도 고양시 덕양구 주교동 600 | 140 | 126.832072 | 37.658042 | POINT (126.83207 37.65804) |
# 주차장 데이터 로드
parking_df = pd.read_csv(os.path.join(DATA_RAW_PATH, "주차장정보.csv"))
# 변수명을 변경한다.
parking_df.columns = ["NAME", "ROAD_NM_ADDR", "LOT_NUM_ADDR", "SPACE", "LON", "LAT"]
parking_df.head()
# LAT, LON 를 geometry point 객체로 변환하고 변수로 추가한다.
import pyproj
from fiona.crs import from_epsg
coulumns = ["NAME", "ROAD_NM_ADDR", "LOT_NUM_ADDR", "SPACE", "LON", "LAT"]
parking_df = gpd.GeoDataFrame(parking_df[coulumns], geometry=gpd.points_from_xy(parking_df.LON, parking_df.LAT))
parking_df.crs = default_crs
parking_df.head()
| code definition | dataframe variable name |
|---|---|
| 도로명주소_건물 정보 | address_building_gdf |
| 도시계획(교통시설) | address_road_gdf |
| 고양시_인도 | sidewalk_gdf |
| BDTYP_CD | BDTYP_MN | BULD_NM1 | BULD_NM2 | GU_CD | GU_NM | DONG_NM | DONG_CD | GRO_FLO_CNT | UND_FLO_CNT | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 01003 | 다가구주택 | None | None | 41281 | 경기도 고양시 덕양구 | 경기도 고양시 덕양구 도내동 | 41281105 | 4 | 0 | MULTIPOLYGON (((126.87256 37.62531, 126.87259 .. | . |
# 건물 데이터 로드
address_building_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '도로명주소_건물.geojson'), encoding='cp949')
# 좌표계 변경
address_building_gdf = address_building_gdf.set_crs(default_crs)
print("address_building_gdf 현재 좌표계 :", address_building_gdf.crs)
# EMD_CD 만들어주기
address_building_gdf["EMD_CD"] = address_building_gdf["SIG_CD"] + address_building_gdf["EMD_CD"]
address_building_gdf.head()
# 변수명 변경하기
address_building_gdf.columns = ["BDTYP_CD", "BULD_NM1", "BULD_NM2", "DONG_CD", "GRO_FLO_CNT", "GU_CD", "UND_FLO_CNT", "geometry"]
address_building_gdf.head()
# 건물 사용 코드 merge
address_building_gdf1 = pd.merge(address_building_gdf, building_usage_code_df, left_on="BDTYP_CD", right_on="code", how="left")
address_building_gdf1 = address_building_gdf1.drop(columns=["code"])
address_building_gdf1.columns = ["BDTYP_CD", "BULD_NM1", "BULD_NM2", "DONG_CD", "GRO_FLO_CNT", "GU_CD", "UND_FLO_CNT", "geometry", "BDTYP_MN"]
address_building_gdf1.head()
# dong_code_df의 code 자료형 변경하기
dong_code_df["code"] = dong_code_df["code"].apply(lambda x: str(x))
# 동 코드 주소 merge
address_building_gdf1 = pd.merge(address_building_gdf1, dong_code_df, left_on="DONG_CD", right_on="code", how="left")
address_building_gdf1 = address_building_gdf1.drop(columns=["code"])
address_building_gdf1.columns = ["BDTYP_CD", "BULD_NM1", "BULD_NM2", "DONG_CD", "GRO_FLO_CNT", "GU_CD", "UND_FLO_CNT", "geometry", "BDTYP_MN", "DONG_NM"]
address_building_gdf1.head()
# gu_code_df code 자료형 변경하기
gu_code_df["code"] = gu_code_df["code"].apply(lambda x: str(x))
# 구 코드 주소 merge
address_building_gdf1 = pd.merge(address_building_gdf1, gu_code_df, left_on="GU_CD", right_on="code", how="left")
address_building_gdf1 = address_building_gdf1.drop(columns=["code"])
address_building_gdf1.columns = ["BDTYP_CD", "BULD_NM1", "BULD_NM2", "DONG_CD", "GRO_FLO_CNT", "GU_CD", "UND_FLO_CNT", "geometry", "BDTYP_MN", "DONG_NM", "GU_NM"]
address_building_gdf1.head()
# 변수 순서 정리
address_building_gdf = address_building_gdf1[["BDTYP_CD", "BDTYP_MN", "BULD_NM1", "BULD_NM2", "GU_CD", "GU_NM", "DONG_NM", "DONG_CD", "GRO_FLO_CNT", "UND_FLO_CNT", "geometry"]]
address_building_gdf.head()
# 간단히 plot해본다.
address_building_gdf.plot()
| 변수 | 의미 |
|---|---|
| RN | 도로명 |
| RDS_MAN_NO | 도로구간 일련번호 |
| RBP_CN | 도로 기점 |
| REP_CN | 도로 종점 |
| BSI_INT | 기초간격 |
| ROAD_BT | 도로 폭 |
| ROAD_LT | 도로 길이 |
| WDR_RD_CD | 광역도로구분코드 |
| ROA_CLS_SE | 도로위계기능구분 |
| RDS_DPN_SE | 도로구간종속구분 |
| RDS_MAN_NO | RN | WDR_RD_CD | ROA_CLS_SE | RDS_DPN_SE | ROAD_BT | ROAD_LT | geometry |
|---|---|---|---|---|---|---|---|
| 21930 | 서울외곽순환고속도로 | 1 | 1 | 0 | 50.0 | 128020.00 | MULTILINESTRING ((126.91306 37.67694, 126.9128... |
# 데이터 로드
address_road_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '도로명주소_도로.geojson'), encoding='cp949')
# 좌표계 변경
address_road_gdf = address_road_gdf.set_crs(default_crs)
print("address_road_gdf 현재 좌표계 :", address_road_gdf.crs)
# 변수명 변경하기
address_road_gdf = address_road_gdf[["RDS_MAN_NO", "RN", "WDR_RD_CD", "ROA_CLS_SE", "RDS_DPN_SE", "ROAD_BT", "ROAD_LT", "geometry"]]
address_road_gdf.head()
# 광역도로 구분 코드(WDR_RD_CD) road_code1_df
road_code1_df
# 도로구간종속구분코드(RDS_DPN_SE) road_code2_df
road_code2_df
# 도로위계구분코드(ROA_CLS_SE) road_code3_df
road_code3_df
# 간단히 plot하기
address_road_gdf.plot(linewidth=1)
| 변수 | 의미 |
|---|---|
| UFID | 인도별 고유코드 |
| 재질 | 인도의 재질을 나타내는 코드 |
| 자전거도로유/무 | 자전거도로의 유무를 나타내는 코드 |
| 종류 | 인도의 종류를 나타내는 코드 |
| SD_CD | BICYCLE | MATERIAL | TYPE | geometry |
|---|---|---|---|---|
0 |TRN0400000001AAEA |BYC001 |SWQ001 |SWK001 |MULTILINESTRING ((126.87388 37.59902, 126.8739...|
sidewalk_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '고양시_인도.geojson'), encoding='cp949')
# 좌표계 변경
sidewalk_gdf = sidewalk_gdf.set_crs(default_crs)
print("sidewalk_gdf 현재 좌표계 :", sidewalk_gdf.crs)
sidewalk_gdf.head()
# 변수명 변경
sidewalk_gdf.columns = ["SD_CD", "MATERIAL", "BICYCLE", "TYPE", "geometry"]
# 변수 순서 변경
sidewalk_gdf = sidewalk_gdf[["SD_CD", "BICYCLE", "MATERIAL", "TYPE", "geometry"]]
sidewalk_gdf.head()
sidewalk_gdf.plot(linewidth=1)
| 파일 | 변수 |
|---|---|
| 고양시덕양구_고도 | duckyang_height_gdf |
| 고양시일산동구_고도 | ilsan_east_height_gdf |
| 고양시일산서구_고도 | ilsan_west_height_gdf |
| 고양시전체_고도 | height_gdf |
# 고도 데이터 로드
duckyang_height_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '고양시_덕양구_고도.geojson'), encoding='cp949')
ilsan_east_height_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '고양시_일산동구_고도.geojson'), encoding='cp949')
ilsan_west_height_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '고양시_일산서구_고도.geojson'), encoding='cp949')
# 좌표계 변경
duckyang_height_gdf = duckyang_height_gdf.set_crs(default_crs)
ilsan_east_height_gdf = ilsan_east_height_gdf.set_crs(default_crs)
ilsan_west_height_gdf = ilsan_west_height_gdf.set_crs(default_crs)
print("duckyang_height_gdf 현재 좌표계 :", duckyang_height_gdf.crs)
print("ilsan_east_height_gdf 현재 좌표계 :", ilsan_east_height_gdf.crs)
print("ilsan_west_height_gdf 현재 좌표계 :", ilsan_west_height_gdf.crs)
# 변수명 변경
duckyang_height_gdf.columns = ["ID", "HEIGHT", "geometry"]
duckyang_height_gdf.head()
# 변수명 변경
ilsan_east_height_gdf.columns = ["ID", "HEIGHT", "geometry"]
ilsan_east_height_gdf.head()
# 변수명 변경
ilsan_west_height_gdf.columns = ["ID", "HEIGHT", "geometry"]
ilsan_west_height_gdf.head()
# 임의의 포인트에 대해 위치 확인
ax = duckyang_height_gdf.plot(color="lightgreen")
ilsan_east_height_gdf.plot(color="black", ax=ax)
ilsan_west_height_gdf.plot(color="skyblue", ax=ax)
plt.show()
# 합해서 plot 해보기
height_gdf = copy.deepcopy(duckyang_height_gdf)
height_gdf = height_gdf.append(ilsan_east_height_gdf)
height_gdf = height_gdf.append(ilsan_west_height_gdf)
height_gdf.reset_index(inplace=True, drop=True)
height_gdf.plot()
| 파일 | 변수 |
|---|---|
| 도시계획(공간시설) | city_plan_space_gdf |
| 도시계획(교통시설) | city_plan_traffic_gdf |
| 용도지역지구(습지보호지역) | limit_area_gdf |
# 도시 계획 데이터 로드
city_plan_space_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '도시계획(공간시설).geojson'), encoding='cp949')
# city_plan_public_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '도시계획(공공문화제육시설).geojson'), encoding='cp949')
city_plan_traffic_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '도시계획(교통시설).geojson'), encoding='cp949')
limit_area_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '용도지역지구(습지보호지역).geojson'), encoding='cp949')
# goyang_geo_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, '고양시 지적도.geojson'), encoding='cp949')
# 좌표계 변경
city_plan_space_gdf = city_plan_space_gdf.set_crs(default_crs)
city_plan_traffic_gdf = city_plan_traffic_gdf.set_crs(default_crs)
limit_area_gdf = limit_area_gdf.set_crs(default_crs)
print("city_plan_space_gdf 현재 좌표계 :", city_plan_space_gdf.crs)
print("city_plan_traffic_gdf 현재 좌표계 :", city_plan_traffic_gdf.crs)
print("limit_area_gdf 현재 좌표계 :", limit_area_gdf.crs)
city_plan_space_gdf.head()
city_plan_space_gdf.plot()
city_plan_traffic_gdf.head()
city_plan_traffic_gdf.plot()
limit_area_gdf.head()
limit_area_gdf.plot()
탐험적 데이터 분석을 통해 해당 문제가 어떠한 특징을 갖고 있고 문제점을 가지고 있는지 이해한다.
# 시각화에 앞서 고양시 중심 위치를 정한다.
center = [37.6559448, 126.8026949]
# 운영이력에서 각 자전거 정류장별 이용량(수요) 추출 함수
def createDemandDf(route_df, unit='all', remove_stations=[0, 998, 999]):
"""
input : route_df , unit = ['all', 'month', 'day']
"""
possible_unit = ['all', 'month', 'day']
if unit not in possible_unit:
print("unit Error!")
return
if type(route_df['LEAS_DATE'][0]) == str:
route_df['LEAS_DATE'] = route_df['LEAS_DATE'].apply(lambda x: datetime.strptime(x, "%Y-%m-%d %H:%M:%S"))
if 'LEAS_DATE_DATE' not in route_df.columns:
route_df['LEAS_DATE_DATE'] = route_df['LEAS_DATE'].apply(lambda x: x.date())
if type(route_df['RTN_DATE'][0]) == str:
route_df['RTN_DATE'] = route_df['RTN_DATE'].apply(lambda x: datetime.strptime(x, "%Y-%m-%d %H:%M:%S"))
if 'RTN_DATE_DATE' not in route_df.columns:
route_df['RTN_DATE_DATE'] = route_df['RTN_DATE'].apply(lambda x: x.date())
now_date = min(route_df['LEAS_DATE'])
end_date = max(route_df['LEAS_DATE'])
stations = sorted(route_df['LEAS_STATION'].unique())
for remove_station in remove_stations:
if remove_station in stations:
stations.remove(remove_station)
stations_num = len(stations)
if unit == 'all':
total_date = (end_date - now_date).days + 1
numpy_route_demand = np.chararray((stations_num, total_date, 2), itemsize=10)
for i in range(len(stations)):
if i == 0:
j = 0
while now_date <= end_date:
numpy_route_demand[i, j] = [str(stations[i]), now_date]
now_date += timedelta(days = 1)
j += 1
else:
temp_numpy_route_demand = copy.deepcopy(numpy_route_demand[0])
temp_numpy_route_demand[:,0] = str(stations[i])
numpy_route_demand[i] = temp_numpy_route_demand
demand_df = pd.DataFrame(columns=['STATION', 'DATE'])
for i in range(len(stations)):
temp_demand_df = pd.DataFrame(numpy_route_demand[i], columns=['STATION', 'DATE'])
demand_df = demand_df.append(temp_demand_df)
demand_df.reset_index(drop=True, inplace=True)
demand_df['STATION'] = demand_df['STATION'].apply(lambda x : int(x))
demand_df['DATE'] = demand_df['DATE'].apply(lambda x : datetime.strptime(x.decode('utf-8'), "%Y-%m-%d"))
demand_df['DATE'] = demand_df['DATE'].apply(lambda x : x.date())
lea_demand_df = route_df[['LEAS_STATION', 'LEAS_DATE_DATE', 'LEAS_NO']].groupby(['LEAS_STATION', 'LEAS_DATE_DATE']).count()
lea_demand_df.reset_index(inplace=True)
lea_demand_df.columns = ['STATION', 'DATE', 'LEAS_COUNTS']
demand_df = pd.merge(demand_df, lea_demand_df, how='left', on=['STATION', 'DATE'])#left_on=['STATION', 'DATE'], right_on=['LEAS_STATION', 'LEAS_DATE_DATE'])
rtn_demand_df = route_df[['RTN_STATION', 'RTN_DATE_DATE', 'LEAS_NO']].groupby(['RTN_STATION', 'RTN_DATE_DATE']).count()
rtn_demand_df.reset_index(inplace=True)
rtn_demand_df.columns = ['STATION', 'DATE', 'RTN_COUNTS']
demand_df = pd.merge(demand_df, rtn_demand_df, how='left', on=['STATION', 'DATE'])
demand_df.fillna(0, inplace=True)
return demand_df
elif unit == 'month':
if 'LEAS_DATE_YEAR' not in route_df.columns:
route_df['LEAS_DATE_YEAR'] = route_df['LEAS_DATE'].apply(lambda x: x.year)
if 'LEAS_DATE_MONTH' not in route_df.columns:
route_df['LEAS_DATE_MONTH'] = route_df['LEAS_DATE'].apply(lambda x: x.month)
if 'RTN_DATE_YEAR' not in route_df.columns:
route_df['RTN_DATE_YEAR'] = route_df['RTN_DATE'].apply(lambda x: x.year)
if 'RTN_DATE_MONTH' not in route_df.columns:
route_df['RTN_DATE_MONTH'] = route_df['RTN_DATE'].apply(lambda x: x.month)
numpy_route_demand = np.chararray((stations_num, 3*12, 3), itemsize=10)
for i in range(len(stations)):
for year in range(2017, 2020):
for month in range(1, 13):
numpy_route_demand[i, (year-2017)*12 + (month-1)] = [str(stations[i]), year, month]
demand_df = pd.DataFrame(columns=['STATION', 'YEAR', 'MONTH'])
for i in range(len(stations)):
temp_demand_df = pd.DataFrame(numpy_route_demand[i], columns=['STATION', 'YEAR', 'MONTH'])
demand_df = demand_df.append(temp_demand_df)
demand_df.reset_index(drop=True, inplace=True)
demand_df['STATION'] = demand_df['STATION'].apply(lambda x : int(x))
demand_df['YEAR'] = demand_df['YEAR'].apply(lambda x : int(x))
demand_df['MONTH'] = demand_df['MONTH'].apply(lambda x : int(x))
lea_demand_df = route_df[['LEAS_STATION', 'LEAS_DATE_YEAR', 'LEAS_DATE_MONTH', 'LEAS_NO']].groupby(['LEAS_STATION', 'LEAS_DATE_YEAR', 'LEAS_DATE_MONTH']).count()
lea_demand_df.reset_index(inplace=True)
lea_demand_df.columns = ['STATION', 'YEAR', 'MONTH', 'LEAS_COUNTS']
demand_df = pd.merge(demand_df, lea_demand_df, how='left', on=['STATION', 'YEAR', 'MONTH'])
rtn_demand_df = route_df[['RTN_STATION', 'RTN_DATE_YEAR', 'RTN_DATE_MONTH', 'LEAS_NO']].groupby(['RTN_STATION', 'RTN_DATE_YEAR', 'RTN_DATE_MONTH']).count()
rtn_demand_df.reset_index(inplace=True)
rtn_demand_df.columns = ['STATION', 'YEAR', 'MONTH', 'RTN_COUNTS']
demand_df = pd.merge(demand_df, rtn_demand_df, how='left', on=['STATION', 'YEAR', 'MONTH'])
demand_df.fillna(0, inplace=True)
return demand_df
elif unit == 'day':
if 'LEAS_DATE_DAY' not in route_df.columns:
route_df['LEAS_DATE_DAY'] = route_df['LEAS_DATE_DATE'].apply(lambda x: x.weekday())
if 'RTN_DATE_DAY' not in route_df.columns:
route_df['RTN_DATE_DAY'] = route_df['RTN_DATE_DATE'].apply(lambda x: x.weekday())
numpy_route_demand = np.chararray((stations_num, 7, 2), itemsize=10)
for i in range(len(stations)):
for day in range(7):
numpy_route_demand[i, day] = [str(stations[i]), day]
demand_df = pd.DataFrame(columns=['STATION', 'DAY'])
for i in range(len(stations)):
temp_demand_df = pd.DataFrame(numpy_route_demand[i], columns=['STATION', 'DAY'])
demand_df = demand_df.append(temp_demand_df)
demand_df.reset_index(drop=True, inplace=True)
demand_df['STATION'] = demand_df['STATION'].apply(lambda x : int(x))
demand_df['DAY'] = demand_df['DAY'].apply(lambda x : int(x))
lea_demand_df = route_df[['LEAS_STATION', 'LEAS_DATE_DAY', 'LEAS_NO']].groupby(['LEAS_STATION', 'LEAS_DATE_DAY']).count()
lea_demand_df.reset_index(inplace=True)
lea_demand_df.columns = ['STATION', 'DAY', 'LEAS_COUNTS']
demand_df = pd.merge(demand_df, lea_demand_df, how='left', on=['STATION', 'DAY'])
rtn_demand_df = route_df[['RTN_STATION', 'RTN_DATE_DAY', 'LEAS_NO']].groupby(['RTN_STATION', 'RTN_DATE_DAY']).count()
rtn_demand_df.reset_index(inplace=True)
rtn_demand_df.columns = ['STATION', 'DAY', 'RTN_COUNTS']
demand_df = pd.merge(demand_df, rtn_demand_df, how='left', on=['STATION', 'DAY'])
demand_df.fillna(0, inplace=True)
return demand_df
import time
def join_height2station(gdf1, gdf2):
"""
input : (gdf_station, gdf_pop_dist)
output : ret_gdf
describe : station gdf에 해당 정류장 기반으로 해발고도를 추가해준다.
"""
ret_df = copy.deepcopy(gdf1)
start = time.time()
for i in range(ret_df.shape[0]):
if ret_df.shape[0] > 1000 & i % 1000 == 0:
print("진행중 : {} / {}".format(i, ret_df.shape[0]))
print("진행시간 : ", time.time()-start)
station = gdf1.iloc[i]["geometry"]
temp_df = gdf2[gdf2.contains(station)]
ret_df.loc[i, "altitude"] = np.mean(list(temp_df["HEIGHT"]))
return ret_df
# 정류장별 이용량 모두 합해주는 함수
def createDemandTotal(x):
df = {}
df['LEAS_COUNTS'] = x['LEAS_COUNTS'].sum()
df['RTN_COUNTS'] = x['RTN_COUNTS'].sum()
df['TOTAL'] = df['LEAS_COUNTS'] + df['RTN_COUNTS']
df['RTN-LEAS'] = df['RTN_COUNTS'] - df['LEAS_COUNTS']
start_index = 0
if 0 in list(x['LEAS_COUNTS']) and 1 in list(x['LEAS_COUNTS']):
start_index = list(np.array(list(x['LEAS_COUNTS'])) > 0).index(True)
df['LEAS_MEAN'] = x['LEAS_COUNTS'][start_index:].mean()
start_index = 0
if 0 in list(x['RTN_COUNTS']) and 1 in list(x['RTN_COUNTS']):
start_index = list(np.array(list(x['RTN_COUNTS'])) > 0).index(True)
df['RTN_MEAN'] = x['RTN_COUNTS'][start_index:].mean()
start_index = 0
if 0 in list(x['TOTAL']) and 1 in list(x['TOTAL']):
start_index = list(np.array(list(x['TOTAL'])) > 0).index(True)
df['TOTAL_MEAN'] = x['TOTAL'][start_index:].mean()
if 0 in list(x['RTN-LEAS']) and 1 in list(x['RTN-LEAS']):
start_index = list(np.array(list(x['RTN-LEAS'])) > 0).index(True)
df['RTN-LEAS_MEAN'] = x['RTN-LEAS'][start_index:].mean()
return pd.Series(df, index=['LEAS_COUNTS', 'RTN_COUNTS', 'TOTAL', 'RTN-LEAS', 'LEAS_MEAN', 'RTN_MEAN', 'TOTAL_MEAN', 'RTN-LEAS_MEAN'])
# 행정동별 자전거 정류장 정보를 모두 합해주는 함수
def calculateStationPerTown(gdf1, gdf2, rel_method=gpd.GeoSeries.within):
"""
input : (gdf_station, gdf_pop_dist, latitude, longitude, coverage, space_relation_method, cal_method)
output : ret_gdf
describe : 행정동 데이터 행 하나에 해당하는 정류장의 개수를 feature로 추가한다.
space_relation_method = gpd.GeoSeries.within, gpd.GeoSeries.intersects
"""
ret_df = copy.deepcopy(gdf1)
for i in range(ret_df.shape[0]):
town = gdf1.iloc[i]["geometry"]
temp_df = gdf2[rel_method(gdf2, town)]
if temp_df.shape[0] == 0:
ret_df.loc[i, "STATION_COUNT"] = 0
ret_df.loc[i, "STATION_RTN_COUNTS"] = 0
ret_df.loc[i, "STATION_LEAS_COUNTS"] = 0
ret_df.loc[i, "STATION_RTN-LEAS"] = 0
ret_df.loc[i, "STATION_TOTAL"] = 0
ret_df.loc[i, "STATION_RTN_MEAN"] = 0
ret_df.loc[i, "STATION_LEAS_MEAN"] = 0
ret_df.loc[i, "STATION_RTN-LEAS_MEAN"] = 0
ret_df.loc[i, "STATION_TOTAL_MEAN"] = 0
else:
ret_df.loc[i, "STATION_COUNT"] = temp_df.shape[0]
ret_df.loc[i, "STATION_RTN_COUNTS"] = np.sum(temp_df["RTN_COUNTS"])
ret_df.loc[i, "STATION_LEAS_COUNTS"] = np.sum(temp_df["LEAS_COUNTS"])
ret_df.loc[i, "STATION_RTN-LEAS"] = np.sum(temp_df["RTN-LEAS"])
ret_df.loc[i, "STATION_TOTAL"] = np.sum(temp_df["TOTAL"])
ret_df.loc[i, "STATION_RTN_MEAN"] = np.mean(temp_df["RTN_MEAN"])
ret_df.loc[i, "STATION_LEAS_MEAN"] = np.mean(temp_df["LEAS_MEAN"])
ret_df.loc[i, "STATION_RTN-LEAS_MEAN"] = np.mean(temp_df["RTN-LEAS_MEAN"])
ret_df.loc[i, "STATION_TOTAL_MEAN"] = np.mean(temp_df["TOTAL_MEAN"])
return ret_df
| STATION_ID | STATION_NM | HOLDER_CAPACITY | LAT | LON | geometry | altitude |
|---|---|---|---|---|---|---|
| 101 | 어울림마을 701동 앞 | 20 | 37.654775 | 126.834584 | POINT (126.83458 37.65477) | 14.0 |
# 사용하기 전 변수 세팅
args = [station_df, height_gdf]
station_df = join_height2station(*args)
station_df.head()
| STATION_ID | STATION_NM | HOLDER_CAPACITY | LAT | LON | geometry | altitude | LEAS_COUNTS | RTN_COUNTS | TOTAL | RTN-LEAS | LEAS_MEAN | RTN_MEAN | TOTAL_MEAN | RTN-LEAS_MEAN |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 101 | 어울림마을 701동 앞 | 20 | 37.654775 | 126.834584 | POINT (941305.151 1961905.875) | 14.0 | 8366.0 | 9207.0 | 17573.0 | 841.0 | 7.640183 |
# 자전거 정류장별 매일 반입 반출량 데이터 만들기
df_station_leasrtn_all = createDemandDf(route_df)
df_station_leasrtn_all['TOTAL'] = df_station_leasrtn_all['LEAS_COUNTS'] + df_station_leasrtn_all['RTN_COUNTS']
df_station_leasrtn_all['RTN-LEAS'] = df_station_leasrtn_all['RTN_COUNTS'] - df_station_leasrtn_all['LEAS_COUNTS']
df_station_leasrtn_all.head()
station_demand_df = df_station_leasrtn_all.groupby(df_station_leasrtn_all['STATION']).apply(createDemandTotal)
station_demand_df = station_demand_df.reset_index()
station_demand_df.head()
# 자전거 거치대 데이터 EDA용 copy해두기
gdf_station = station_df.copy(deep=True)
gdf_station = gdf_station.to_crs(meter_crs) # 좌표계 변환
# 100m*100m 인구분포 데이터 EDA용 copy해두기
gdf_pop_dist = population_dist_gdf.copy(deep=True)
# 좌표계 설정
gdf_pop_dist.crs = default_crs
gdf_pop_dist = gdf_pop_dist.to_crs(meter_crs)
print("gdf_pop_dist 현재 좌표계 :", gdf_pop_dist.crs)
# na 처리
gdf_pop_dist["val"] = gdf_pop_dist["val"].fillna(0)
gdf_pop_dist.head()
운행이력에 있는 데이터에 있는 총 정류장 갯수를 조사한 결과 155개이고,
정류장에 대한 정보를 불러왔을 때의 정류장 갯수는 163개이다.
운행이력에 담겨있지 않은 데이터는 통계치를 산정하는데 문제가 발생하기 때문에 정류소에 대한 정보를 취합할 때 제거한다.
# 자전거 정류장 데이터에 정류장에 대한 통계치 붙이기(merge)
gdf_station = pd.merge(gdf_station, station_demand_df, left_on='STATION_ID', right_on="STATION")
del gdf_station["STATION"]
gdf_station.head()
| GU_CD | GU_NM | DONG_CD | DONG_NM | geometry | STATION_COUNT | STATION_RTN_COUNTS | STATION_LEAS_COUNTS | STATION_RTN-LEAS | STATION_TOTAL | STATION_RTN_MEAN | STATION_LEAS_MEAN | STATION_RTN-LEAS_MEAN | STATION_TOTAL_MEAN |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 31101 | 덕양구 | 3110151 | 주교동 | MULTIPOLYGON (((941371.354 1964353.054, 941373... | 1.0 | 2950.0 | 3671.0 | -721.0 6621.0 | 2.696527 | 3.355576 | -0.659049 | 6.052102 |
# 행정동 경계 데이터 EDA용도로 복사한다.
gdf_boundary_town = administration_boundary_town_gdf.copy(deep=True)
# 좌표계 변환
gdf_boundary_town = gdf_boundary_town.to_crs(meter_crs)
print("변환 좌표계 :", gdf_boundary_town.crs)
gdf_boundary_town.head()
# 사용하기 전 변수 세팅
args = [gdf_boundary_town, gdf_station, gpd.GeoSeries.within]
# 함수 실행
gdf_boundary_town = calculateStationPerTown(*args)
gdf_boundary_town.head()
# 조회기준 2019년 12월 데이터만 사용
# EMD_KOR_NM에서 행정동만 골라내기
gdf_town_pop = population_statistics_df.copy(deep=True)
gdf_town_pop = gdf_town_pop[gdf_town_pop["CHECK_DATE"] == "2019년12월"].reset_index(drop=True)
gdf_town_pop = gdf_town_pop.drop(index=[0,1,22,34]).reset_index(drop=True)
gdf_town_pop = gdf_town_pop.drop('CHECK_DATE', axis=1)
gdf_town_pop["DONG_NM"] = gdf_town_pop["EMD_KOR_NM"].apply(lambda x: x.split(' ')[3].split('(')[0])
del gdf_town_pop["EMD_KOR_NM"]
# 행정동별 정거장 데이터와 인구통계 데이터 합치기
gdf_boundary_town = pd.merge(gdf_boundary_town, gdf_town_pop, on = "DONG_NM", how = "left")
# 정거장당 인구 수 변수 추가
gdf_boundary_town["STATION_PER_POP"] = gdf_boundary_town["TOTAL_POP"] / gdf_boundary_town["STATION_COUNT"]
gdf_boundary_town["RIDE_PER_POP"] = gdf_boundary_town["STATION_TOTAL"] / gdf_boundary_town["TOTAL_POP"]
# 정거장이 0인 경우, 0으로 처리
for i in range(len(gdf_boundary_town)):
if np.isinf(gdf_boundary_town["STATION_PER_POP"][i]):
gdf_boundary_town["STATION_PER_POP"][i] = 0
gdf_boundary_town.head()
# folium 마커의 모양을 html로 설정한다.
table = """
<!DOCTYPE html>
<html>
<head>
<style>
table {{
width:100%;
}}
table, th, td {{
border: 1px solid black;
border-collapse: collapse;
}}
th, td {{
padding: 5px;
text-align: left;
}}
table#t01 tr:nth-child(odd) {{
background-color: #eee;
}}
table#t01 tr:nth-child(even) {{
background-color:#fff;
}}
</style>
</head>
<body>
<table id="t01">
<tr>
<td>Name</td>
<td>{}</td>
</tr>
</table>
</body>
</html>
""".format
# folium rendering 관계로 버스정류장은 100개만 plot한다.
# colab에서는 모두 렌더링이되는 것을 확인했으나, 구동환경의 차이로 불가한 것을 확인했다.
import folium
from folium import IFrame
from folium.plugins import MarkerCluster, FeatureGroupSubGroup
# 행사장 공간정보 eventhall_space_df
# 공연장 공간정보 concerthall_space_df
# 박물관/미술관 공간정보 museum_space_df
# 버스정류장 정보 busstop_df
# 지하철 역 정보 subway_space_df
# 체육시설 정보 sport_complex_info_df
# 주차장 정보 parking_df
m = folium.Map(location=center, tiles="OpenStreetMap", zoom_start=12)
fg = folium.FeatureGroup(name='ALL')
m.add_child(fg)
g1 = FeatureGroupSubGroup(fg, 'bike station[star]')
m.add_child(g1)
g2 = FeatureGroupSubGroup(fg, 'event hall[cloud]')
m.add_child(g2)
g3 = FeatureGroupSubGroup(fg, 'concert hall[heart]')
m.add_child(g3)
g4 = FeatureGroupSubGroup(fg, 'museum[lock]')
m.add_child(g4)
g5 = FeatureGroupSubGroup(fg, 'busstop[ok]')
m.add_child(g5)
g6 = FeatureGroupSubGroup(fg, 'subway station[tag]')
m.add_child(g6)
g7 = FeatureGroupSubGroup(fg, 'sport complex[tint]')
m.add_child(g7)
g8 = FeatureGroupSubGroup(fg, 'parking lot[camera]')
m.add_child(g8)
width, height = 310,150
popups, locations = [], []
marker_cluster = MarkerCluster().add_to(g1)
for idx, row in station_df.iterrows():
folium.Marker(
location=[row['geometry'].y, row['geometry'].x],
popup=IFrame(table(row['STATION_NM']), width=width, height=height),
icon=folium.Icon(color='red',icon='star'),
).add_to(marker_cluster)
popups, locations = [], []
marker_cluster = MarkerCluster().add_to(g2)
for idx, row in eventhall_space_df.iterrows():
folium.Marker(
location=[row['geometry'].y, row['geometry'].x],
popup=IFrame(table(row['NAME']), width=width, height=height),
icon=folium.Icon(color='blue',icon='cloud'),
).add_to(marker_cluster)
popups, locations = [], []
marker_cluster = MarkerCluster().add_to(g3)
for idx, row in concerthall_space_df.iterrows():
folium.Marker(
location=[row['geometry'].y, row['geometry'].x],
popup=IFrame(table(row['NAME']), width=width, height=height),
icon=folium.Icon(color='purple',icon='heart'),
).add_to(marker_cluster)
popups, locations = [], []
marker_cluster = MarkerCluster().add_to(g4)
for idx, row in museum_space_df.iterrows():
folium.Marker(
location=[row['geometry'].y, row['geometry'].x],
popup=IFrame(table(row['NAME']), width=width, height=height),
icon=folium.Icon(color='magenda',icon='lock'),
).add_to(marker_cluster)
popups, locations = [], []
marker_cluster = MarkerCluster().add_to(g5)
for idx, row in busstop_df.iterrows():
if idx == 100: break
folium.Marker(
location=[row['geometry'].y, row['geometry'].x],
popup=IFrame(table(row['STATION_NM']), width=width, height=height),
icon=folium.Icon(color='orange',icon='ok'),
).add_to(marker_cluster)
popups, locations = [], []
marker_cluster = MarkerCluster().add_to(g6)
for idx, row in subway_space_df.iterrows():
folium.Marker(
location=[row['geometry'].y, row['geometry'].x],
popup=IFrame(table(row['STATION_NM']), width=width, height=height),
icon=folium.Icon(color='green',icon='tag'),
).add_to(marker_cluster)
popups, locations = [], []
marker_cluster = MarkerCluster().add_to(g7)
for idx, row in sport_complex_info_df.iterrows():
folium.Marker(
location=[row['geometry'].y, row['geometry'].x],
popup=IFrame(table(row['NAME']), width=width, height=height),
icon=folium.Icon(color='yellow',icon='tint'),
).add_to(marker_cluster)
popups, locations = [], []
marker_cluster = MarkerCluster().add_to(g8)
for idx, row in parking_df.iterrows():
folium.Marker(
location=[row['geometry'].y, row['geometry'].x],
popup=IFrame(table(row['NAME']), width=width, height=height),
icon=folium.Icon(color='red',icon='camera'),
).add_to(marker_cluster)
folium.GeoJson(
administration_boundary_town_gdf,
name='Administration Boundary Layer',
style_function=lambda feature: {
'fillOpacity': 0.2,
'weight': 2,
"color": "#1C1C1C",
'lineColor': '#E2E2E2',
"opacity": 0.4,
}
).add_to(m)
folium.LayerControl(collapsed=False).add_to(m)
m
대부분의 자전거 정류장은 도시지역에 분포한다는 사실을 알 수 있다.
# 행정동 인구수를 파악한다.
from branca.colormap import linear
m = folium.Map(center, zoom_start=12,tiles=None)
# feature groups
# feature_group0 = folium.FeatureGroup(name='TOTAL_POP',overlay=False).add_to(m) # 행정동 인구수
# feature_group1= folium.FeatureGroup(name='STATION_COUNT',overlay=False).add_to(m) # 행정동 내 자전거 정류장 수
feature_group2 = folium.FeatureGroup(name='STATION_PER_POP',overlay=False).add_to(m) # 자전거 정류장 하나 당 커버하는 행정동 인구수
fs = [feature_group2]#,feature_group1]#,feature_group2]
feature_list = ["STATION_PER_POP"]#,"STATION_COUNT"]#,"STATION_PER_POP"]
for i in range(len(feature_list)):
choropleth1 = folium.Choropleth(
geo_data=gdf_boundary_town,
name='choropleth',
data=gdf_boundary_town,
columns=['DONG_NM', feature_list[i]],
key_on='feature.properties.DONG_NM',
fill_color='YlGn',
nan_fill_color="black",
fill_opacity=0.7,
line_opacity=0.2,
legend_name='Arrival (in Quintal)',
highlight=True,
line_color='black').geojson.add_to(fs[i])
colormap = linear.YlGn_09.scale(
gdf_boundary_town[feature_list[i]].min(),
gdf_boundary_town[feature_list[i]].max()).to_step(10)
colormap.caption = feature_list[i]
colormap.add_to(m)
folium.TileLayer('cartodbpositron',overlay=True,name="background").add_to(m)
folium.LayerControl('topright',collapsed=False).add_to(m)
m
# 행정동 내 인구 대비 자전거 거치대 수를 파악한다.
from branca.colormap import linear
m = folium.Map(center, zoom_start=12,tiles=None)
# feature groups
# feature_group0 = folium.FeatureGroup(name='TOTAL_POP',overlay=False).add_to(m) # 행정동 인구수
# feature_group1= folium.FeatureGroup(name='STATION_COUNT',overlay=False).add_to(m) # 행정동 내 자전거 정류장 수
feature_group2 = folium.FeatureGroup(name='STATION_PER_POP',overlay=False).add_to(m) # 자전거 정류장 하나 당 커버하는 행정동 인구수
fs = [feature_group2]#,feature_group1]#,feature_group2]
feature_list = ["STATION_PER_POP"]#,"STATION_COUNT"]#,"STATION_PER_POP"]
for i in range(len(feature_list)):
choropleth1 = folium.Choropleth(
geo_data=gdf_boundary_town,
name='choropleth',
data=gdf_boundary_town,
columns=['DONG_NM', feature_list[i]],
key_on='feature.properties.DONG_NM',
fill_color='YlGn',
nan_fill_color="black",
fill_opacity=0.7,
line_opacity=0.2,
legend_name='Arrival (in Quintal)',
highlight=True,
line_color='black').geojson.add_to(fs[i])
colormap = linear.YlGn_09.scale(
gdf_boundary_town[feature_list[i]].min(),
gdf_boundary_town[feature_list[i]].max()).to_step(10)
colormap.caption = feature_list[i]
colormap.add_to(m)
folium.TileLayer('cartodbpositron',overlay=True,name="background").add_to(m)
folium.LayerControl('topright',collapsed=False).add_to(m)
m
인구 수가 많은 주엽동, 마두동의 경우 충분한 자전거 정류장이 존재했기에 상대적으로 낮은 거치대 대비 인구 수를 보이나,
행신동은 인구 대비 거치대 수가 적고, 고양시 북부의 경우 아예 거치대가 존재하지 않는 등 배치의 불균형성이 보인다.
# 100m*100m 인구분포 데이터 전처리
# 좌표계 변환
population_dist_gdfM = population_dist_gdf.copy(deep=True)
population_dist_gdfM = population_dist_gdfM.to_crs(meter_crs)
print("변환 좌표계 :", population_dist_gdfM.crs)
# Nan = 0으로 변환
population_dist_gdfM["val"] = population_dist_gdfM["val"].fillna(0)
# val에 따라 데이터 분리
population_dist_gdfM_not_zero = population_dist_gdfM[population_dist_gdfM["val"] > 0]
population_dist_gdfM_zero = population_dist_gdfM[population_dist_gdfM["val"] == 0]
population_dist_gdfM_not_zero = population_dist_gdfM_not_zero.reset_index(drop=True)
population_dist_gdfM_not_zero.head()
# 인구분포도 바탕으로 자전거 거치대의 200m 반경 커버 지역 시각화
station_boundary_200 = gpd.GeoDataFrame({'geometry': gdf_station.buffer(200)})
origin_ = population_dist_gdfM_not_zero.groupby(['gid']).apply(lambda gr : gr.area.sum())
ax = gdf_boundary_town.plot(color="silver", figsize=(16,16))
population_dist_gdfM_not_zero.plot(column='val', legend=True, scheme='Fisher_Jenks', cmap="copper", k=10, ax=ax)
gdf_station.plot(ax=ax, marker='v', color='black')
station_boundary_200.boundary.plot(ax=ax,linewidth = 0.01, facecolor = 'white', color='red')
ax.set_axis_off()
ax.set_title("bike station coverage based on population distribution")
plt.show()
시각화 결과, 거주 인구가 상당히 있음에도 불구하고 거치대가 설치되지 않은 곳이 많은 것을 알 수 있다. 우리는 이 부분에서 미배치된 지역과, 배치된 지역을 분리하여 추후에 사용하기로 하였다.
# 100x100 인구분포도를 기반으로 자전거 거치대의 200m 반경 미커버 지역 비율 계산
# 현재 거치대가 cover하지 못하는 gid의 비율을 계산한다.
dif_area = gpd.overlay(population_dist_gdfM_not_zero, station_boundary_200, how='difference')
dif_area = dif_area.dissolve(by='gid')
# 커버 지역의 총 인구수, 인구 비율 (인구분포도 기반)
# 100% cover하는 지역의 return값이 nan으로 생성된다. 이를 0으로 변환한다.
cover_gid_p = round(dif_area.area / origin_ * 100)
cover_gid_p = cover_gid_p.fillna(0)
cover_gid_p[cover_gid_p.isna()]
cover_gid_p.head()
# 미커비비율 데이터 만들기
cover_gid_p = cover_gid_p.reset_index()
cover_gid_p.columns = ["gid", "not_cover_percentage"]
cover_gid_p.head()
# 미커버 비율에 대한 테이블을 원래 100x100 인구분포도 데이터와 gid를 기준으로 merge한다.
# val이 0보다 큰 인구분포 데이터(거주 인구가 0이 아닌 지역)와 합치기(merge)
population_dist_gdfM_not_zero_coverage = population_dist_gdfM_not_zero.copy(deep=True)
population_dist_gdfM_not_zero_coverage["geometry"] = population_dist_gdfM_not_zero_coverage["geometry"].to_crs(meter_crs)
population_dist_gdfM_not_zero_coverage = pd.merge(population_dist_gdfM_not_zero, cover_gid_p, on="gid")
population_dist_gdfM_not_zero_coverage.head()
# 미커버비율이 100%인 gid만 골내서 val(인구 수)로 sorting한다.
# 이 부분은 추후에 후보군 선정에서 사용한다.
population_dist_gdfM_not_zero_coverage = population_dist_gdfM_not_zero_coverage[population_dist_gdfM_not_zero_coverage["not_cover_percentage"] == 100]
population_dist_gdfM_not_zero_coverage.sort_values(by=['val'], axis=0, ascending=False, inplace=True)
population_dist_gdfM_not_zero_coverage = population_dist_gdfM_not_zero_coverage.reset_index(drop=True)
population_dist_gdfM_not_zero_coverage.head()
현재 거치대는 두가지 위치 특징을 기반으로 설치되었다고 가정한다.
이러한 가정을 설정한다면, 현재 거치대가 설치된 주변 건물의 특징을 조사한다면, 거치대가 상업지구에 설치되었는지, 거주지구에 설치되었는지 알 수 있다.
우리는 이 단계에서 상업지구와 거주지구를 나누기 위한 변수 "house_counts_mean"를 사용할 것이다. "house_counts_mean"는 거치대 주변 200m내에 존재하는 거주 건물의 수의 평균이다. 만약 200m 반경 내에 4개의 100x100m 격자가 들어가고, 이 격자 각각에 2개 씩의 거주 건물이 들어가 있다면 (4x2)/4 = 2로 생성된다. 이를 적용했을 때, 1을 기준으로 상업지구와 거주인구가 나뉘는 것을 확인할 수 있었다. 이를 좀더 명시적으로 파악하기 위해 네이버 지도에서 해당 위치와의 비교를 통해 검증하였다.

과정은 다음과 같다.
# 100m*100m 빌딩연면적 데이터 전처리
# 좌표계 변환
building_area_gdfM = building_area_gdf.copy(deep=True)
building_area_gdfM = building_area_gdfM.to_crs(meter_crs)
print("변환 좌표계 :", building_area_gdfM.crs)
# Nan = 0으로 변환
building_area_gdfM["val"] = building_area_gdfM["val"].fillna(0)
# val에 따라 데이터 분리
building_area_gdfM_not_zero = building_area_gdfM[building_area_gdfM["val"] > 0]
building_area_gdfM_zero = building_area_gdfM[building_area_gdfM["val"] == 0]
building_area_gdfM_not_zero = building_area_gdfM_not_zero.reset_index(drop=True)
building_area_gdfM_not_zero.head()
# 100m*100m 인구분포 데이터와 빌딩연면적 데이터 합치기(merge)
building_pop_gdfM = pd.merge(building_area_gdfM, population_dist_gdfM, on = "gid", how = "left")
building_pop_gdfM = building_pop_gdfM.drop(columns=["geometry_y"])
building_pop_gdfM.columns = ["gid", "b_area", "geometry", "pop"]
building_pop_gdfM[245:251]
# 인구 수와 빌딩 면적이 0보다 큰 데이터셋으로 분류
building_pop_not_all_zero_gdf = building_pop_gdfM[(building_pop_gdfM["b_area"] > 0) & (building_pop_gdfM["pop"] > 0)]
building_pop_not_zero_gdf = building_pop_gdfM[(building_pop_gdfM["b_area"] > 0) | (building_pop_gdfM["pop"] > 0)]
print(building_pop_not_zero_gdf.shape)
print(building_pop_not_zero_gdf[building_pop_not_zero_gdf["b_area"] == 0].shape)
print(building_pop_not_zero_gdf[building_pop_not_zero_gdf["pop"] == 0].shape)
print(building_pop_not_all_zero_gdf.shape)
building_pop_not_zero_gdf.head()
도로명주소의 건축물용도코드를 활용하여 주거지역을 판단한다.
주거지 관련 건축물용도코드: "01000","01001","01002","01003","01004","02000","02001","02002","02003","02004","02005","02006","02007"
address_building_gdf_temp = address_building_gdf[address_building_gdf["BDTYP_CD"].isin(["01000","01001","01002","01003","01004","02000","02001","02002","02003","02004","02005","02006","02007"])]
address_building_gdf_temp = address_building_gdf_temp.reset_index(drop=True)
address_building_gdf_temp.head()
# # 앞에서 만든 100m*100m 인구분포, 빌딩 면적 데이터의 각각 gid 안에 몇개의 주거용도 건축물이 있는지 count
# import time
# start = time.time()
# num = len(address_building_gdf_temp)
# address_building_gdf_temp['HOUSE_COUNTS'] = 0
# for idx, row in address_building_gdf_temp.iterrows():
# if idx % 1000 == 0:
# print("진행률 : {} / {}".format(idx, num))
# print("진행시간 : {}".format(time.time() - start))
# address_building_gdf_temp['HOUSE_COUNTS'] += gpd.GeoSeries.contains(building_pop_not_zero_gdf['geometry'], row.loc['geometry'].centroid)
house_counts_gdf = gpd.read_file(os.path.join(DATA_RAW_PATH, 'temp.geojson'), encoding='cp949')
# house_counts_gdf = address_building_gdf_temp.copy(deep=True)
house_counts_gdf.head()
# house_counts의 모든 경우를 살펴본다.
house_counts_gdf["HOUSE_COUNTS"].unique()
house_counts_gdf.head()
# 각 gid(격자)에 주거용도 건축물의 수 붙여주기
building_pop_not_zero_gdf= building_pop_not_zero_gdf.reset_index(drop=True)
house_counts_gdf = pd.merge(building_pop_not_zero_gdf, house_counts_gdf, on = "gid", how="left")
# 변수 처리
del house_counts_gdf["geometry_y"]
del house_counts_gdf["b_area_y"]
del house_counts_gdf["pop_y"]
del house_counts_gdf["id"]
house_counts_gdf.columns = ["gid", "b_area", "geometry", "pop", "HOUSE_COUNTS"]
house_counts_gdf.head()
# gid 별 현재 자전거 정류장이 미커버하는 비율 붙여주기
house_counts_cover_gdf = pd.merge(cover_gid_p, house_counts_gdf, on = "gid", how="left")
house_counts_cover_gdf = gpd.GeoDataFrame(house_counts_cover_gdf, geometry="geometry")
# 좌표계 변환
house_counts_cover_gdf["geometry"] = house_counts_cover_gdf["geometry"].to_crs(meter_crs)
print(house_counts_cover_gdf.crs)
house_counts_cover_gdf.head()
# 주거용도 건축물이 1 이상 있는 데이터셋
house_counts_not_zero_gdf = house_counts_gdf[house_counts_gdf["HOUSE_COUNTS"] > 0]
import time
def join_info2station(gdf1, gdf2, coverage, rel_method=gpd.GeoSeries.within, cal_method=sum):
"""
input : (gdf_station, house_counts_cover_gdf, latitude, longitude, coverage, space_relation_method, cal_method)
output : ret_gdf
describe : station gdf에 해당 정류장 기반으로 coverage 안에 들어오는 인구 분포, 면적 분포, 주거용도 건축물의 수를 feature로 추가한 gdf를 반환한다.
space_relation_method = gpd.GeoSeries.within, gpd.GeoSeries.intersects
cal_method = sum, np.mean
"""
ret_df = copy.deepcopy(gdf1)
start = time.time()
for i in range(ret_df.shape[0]):
if ret_df.shape[0] > 1000 and i % 1000 == 0:
print("진행중 {} / {}".format(i, ret_df.shape[0]))
print("진행 시간", time.time()-start)
station = gdf1.iloc[i]["geometry"]
station = station.buffer(coverage)
temp_df = gdf2[rel_method(gdf2, station)]
ret_df.loc[i, f"pop_{coverage}_{cal_method.__name__}"] = cal_method(temp_df["pop"])
ret_df.loc[i, f"b_area_{coverage}_{cal_method.__name__}"] = cal_method(temp_df["b_area"])
ret_df.loc[i, f"HOUSE_COUNTS_{coverage}_mean"] = np.mean(temp_df["HOUSE_COUNTS"])
return ret_df
# 사용하기 전 변수 세팅
args = [gdf_station, house_counts_cover_gdf, 200, gpd.GeoSeries.within, np.sum]
# 함수 실행
gdf_station_cover = join_info2station(*args)
gdf_station_cover.head() # 결과 확인
sns.distplot(gdf_station_cover["HOUSE_COUNTS_200_mean"])
200m 단위로 house_counts_mean을 보았을 때, 0과 1에서 가장 많은 양상을 보였다.
우리는 이 점을 바탕으로 1을 기준으로 상업지구와 주거지구를 나누고자 했다.
# 정류장 주변 주거용 건축물이 평균적으로 1 이하/초과인 데이터셋을 추출하여 확인한다.
station_cover_below1 = gdf_station_cover[gdf_station_cover["HOUSE_COUNTS_200_mean"] <= 1]
station_cover_over1 = gdf_station_cover[gdf_station_cover["HOUSE_COUNTS_200_mean"] > 1]
print(station_cover_below1.shape)
print(station_cover_over1.shape)
station_cover_below1.head()
# 시각화한 인구분포도에서 상업지역 걸러내기
station_boundary_200 = gpd.GeoDataFrame({'geometry': station_cover_below1.buffer(200)})
origin_ = house_counts_gdf.groupby(['gid']).apply(lambda gr : gr.area.sum())
dif_area = gpd.overlay(house_counts_gdf, station_boundary_200, how='difference')
dif_area = dif_area.dissolve(by='gid')
# 상업지역 gid 골라내기
commercial_area_gid_p = round(dif_area.area / origin_ * 100)
commercial_area_gid_p = commercial_area_gid_p.fillna(0)
# 미커비비율 데이터 만들기
commercial_area_gid_p = commercial_area_gid_p.reset_index()
commercial_area_gid_p.columns = ["gid", "not_commerical_area_percentage"]
print(commercial_area_gid_p.shape)
# house counts 데이터에 합치기
house_counts_gdf = pd.merge(house_counts_gdf, commercial_area_gid_p, on="gid", how="left")
print(house_counts_gdf.shape)
house_counts_gdf.head()
# 시각화한 인구분포도에서 주거지역 걸러내기
station_boundary_200 = gpd.GeoDataFrame({'geometry': station_cover_over1.buffer(200)})
origin_ = house_counts_gdf.groupby(['gid']).apply(lambda gr : gr.area.sum())
dif_area = gpd.overlay(house_counts_gdf, station_boundary_200, how='difference')
dif_area = dif_area.dissolve(by='gid')
# 주거지역 gid 골라내기
residential_area_gid_p = round(dif_area.area / origin_ * 100)
residential_area_gid_p = residential_area_gid_p.fillna(0)
# 미커비비율 데이터 만들기
residential_area_gid_p = residential_area_gid_p.reset_index()
residential_area_gid_p.columns = ["gid", "not_residential_area_percentage"]
print(residential_area_gid_p.shape)
# house counts 데이터에 합치기
house_counts_gdf = pd.merge(house_counts_gdf, residential_area_gid_p, on="gid", how="left")
print(house_counts_gdf.shape)
house_counts_gdf.head()
# 좌표계 변경해주기
building_area_zero_gdf["geometry"] = building_area_zero_gdf["geometry"].to_crs(meter_crs)
# 선별된 상업지역 시각화
ax = house_counts_not_zero_gdf.plot( column="pop", cmap='copper', figsize = (30,30), legend=True)
building_area_zero_gdf.plot(color="silver", ax=ax)
station_cover_below1.plot(ax=ax, marker='v', markersize = 30, color='blue')
house_counts_gdf[house_counts_gdf["not_commerical_area_percentage"] < 100].plot(marker='v', color = 'yellow', ax=ax)
ax.set_title("House counts Distribution", fontsize=40)
ax.set_axis_off()
plt.show()
흰 지역은 주거용 건축물이 하나도 없는 지역을 의미한다.
빨간 점으로 표시된 주변 지역은 주거용 건축물이 거의 없으면, 실제 지도와 비교 분석해본 결과 고양시 주요 상업지구로 판단해도 문제가 없다는 것을 확인했다.

# 선별된 주거지역 시각화
ax = house_counts_not_zero_gdf.plot( column="pop", cmap='copper', figsize = (30,30), legend=True)
building_area_zero_gdf.plot(color="silver", ax=ax)
house_counts_gdf[house_counts_gdf["not_residential_area_percentage"] < 100].plot(marker='v', color = 'yellow', ax=ax)
ax.set_title("House counts Distribution", fontsize=20)
ax.set_axis_off()
plt.show()
house_counts_gdf["commerical_area"] = house_counts_gdf["not_commerical_area_percentage"].apply(lambda x: 0 if x == 100.0 else 1 )
house_counts_gdf["residential_area"] = house_counts_gdf["not_residential_area_percentage"].apply(lambda x: 0 if x == 100.0 else 2)
house_counts_gdf["areatype"] = house_counts_gdf["commerical_area"] + house_counts_gdf["residential_area"]
house_counts_gdf.head()
주거 지역 역시 파악해 본 결과 대다수의 거치대가 아파트 근처이거나 거주지임을 확인했다.
sns.set(style="white", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=1,ncols=3)
fig.set_size_inches(12, 10)
sns.boxplot(x="areatype", y="HOUSE_COUNTS", data=house_counts_gdf, ax=axes[0])
sns.boxplot(x="areatype", y="pop", data=house_counts_gdf, ax=axes[1])
sns.boxplot(x="areatype", y="b_area", data=house_counts_gdf, ax=axes[2])
plt.tight_layout()
plt.show()
| areatype | describe |
|---|---|
| 0 | 어느곳에도 포함되지 않는 지역 |
| 1 | 상업지역 |
| 2 | 주거지역 |
| 3 | 상업지역 & 주거지역 |
위에 그래프들을 보면 주거지역이나 상업지역 어느 곳에도 들지 못한 areatype = 0인 구역에 상당히 많은 아웃라이어들이 보인다.
이는 현재 배치된 자전거 거치대을 바탕으로 주거, 상업지역을 판단했기 때문으로, 현재 자전거 정류장이 cover하지 못한 주거 지역들이 다수 존재하는 것으로 판단할 수 있다.
따라서 상업지역은 그 지역의 특성을 대변하는 것으로 사용해도 되지만, 주거지역은 이 데이터만을 믿고 추출하는 것에 무리가 있다. 따라서 추후 주거지역을 기반으로 특성을 추출할 때, 이 데이터 이외에 다른 지표를 기반으로 주거지역을 대변할 필요가 있다.
자전거의 이용횟수를 수요라고 보았을 때, 다양한 관점으로 이를 관찰하며 어떠한 특징이 있는지 확인한다.
# 행정동 내 자전거 거치대 개수를 파악한다.
from branca.colormap import linear
m = folium.Map(center, zoom_start=12,tiles=None)
# feature groups
# feature_group0 = folium.FeatureGroup(name='TOTAL_POP',overlay=False).add_to(m) # 행정동 인구수
feature_group1= folium.FeatureGroup(name='STATION_COUNT',overlay=False).add_to(m) # 행정동 내 자전거 정류장 수
# feature_group2 = folium.FeatureGroup(name='STATION_PER_POP',overlay=False).add_to(m) # 자전거 정류장 하나 당 커버하는 행정동 인구수
fs = [feature_group1]#,feature_group1]#,feature_group2]
feature_list = ["STATION_COUNT"]#,"STATION_COUNT"]#,"STATION_PER_POP"]
for i in range(len(feature_list)):
choropleth1 = folium.Choropleth(
geo_data=gdf_boundary_town,
name='choropleth',
data=gdf_boundary_town,
columns=['DONG_NM', feature_list[i]],
key_on='feature.properties.DONG_NM',
fill_color='YlGn',
nan_fill_color="black",
fill_opacity=0.7,
line_opacity=0.2,
legend_name='Arrival (in Quintal)',
highlight=True,
line_color='black').geojson.add_to(fs[i])
colormap = linear.YlGn_09.scale(
gdf_boundary_town[feature_list[i]].min(),
gdf_boundary_town[feature_list[i]].max()).to_step(10)
colormap.caption = feature_list[i]
colormap.add_to(m)
folium.TileLayer('cartodbpositron',overlay=True,name="background").add_to(m)
folium.LayerControl('topright',collapsed=False).add_to(m)
m
# 행정동별 총 자전거 이용량 (반입량+반출량)의 일별 평균 : 절대 수요을 시각화한다.
from branca.colormap import linear
m = folium.Map(center, zoom_start=12,tiles=None)
# feature groups
feature_group0 = folium.FeatureGroup(name='STATION_TOTAL_MEAN',overlay=False).add_to(m) # 행정동별 총 자전거 이용량
# feature_group1= folium.FeatureGroup(name='STATION_RTN-LEAS_MEAN',overlay=False).add_to(m) # 행정동별 정류장들의 (자전거 반입량 - 반출량) 평균
# feature_group2 = folium.FeatureGroup(name='STATION_RTN_COUNTS',overlay=False).add_to(m) # 행정동별 총 자전거 반입량
# feature_group3 = folium.FeatureGroup(name='STATION_LEAS_COUNTS',overlay=False).add_to(m) # 행정동별 총 자전거 반출량
# feature_group4 = folium.FeatureGroup(name='RIDE_PER_POP',overlay=False).add_to(m) # 행정동별 한명당 자전거 이용량
fs = [feature_group0]#,feature_group1,feature_group2,feature_group3,feature_group4]
feature_list = ["STATION_TOTAL_MEAN"]#,"STATION_RTN-LEAS_MEAN","STATION_RTN_COUNTS","STATION_LEAS_COUNTS","RIDE_PER_POP"]
for i in range(len(feature_list)):
choropleth1 = folium.Choropleth(
geo_data=gdf_boundary_town,
name='choropleth',
data=gdf_boundary_town,
columns=['DONG_NM', feature_list[i]],
key_on='feature.properties.DONG_NM',
fill_color='YlGn',
nan_fill_color="black",
fill_opacity=0.7,
line_opacity=0.2,
legend_name='Arrival (in Quintal)',
highlight=True,
line_color='black').geojson.add_to(fs[i])
colormap = linear.YlGn_09.scale(
gdf_boundary_town[feature_list[i]].min(),
gdf_boundary_town[feature_list[i]].max()).to_step(10)
colormap.caption = feature_list[i]
colormap.add_to(m)
folium.TileLayer('cartodbpositron',overlay=True,name="background").add_to(m)
folium.LayerControl('topright',collapsed=False).add_to(m)
m
# 행정동별 정류장들의 총 자전거 이용량 (반입량-반출량)의 일별 평균 : 상대 수요을 시각화한다.
from branca.colormap import linear
m = folium.Map(center, zoom_start=12,tiles=None)
# feature groups
# feature_group0 = folium.FeatureGroup(name='STATION_TOTAL',overlay=False).add_to(m) # 행정동별 총 자전거 이용량
feature_group1= folium.FeatureGroup(name='STATION_RTN-LEAS_MEAN',overlay=False).add_to(m) # 행정동별 정류장들의 (자전거 반입량 - 반출량) 평균
# feature_group2 = folium.FeatureGroup(name='STATION_RTN_COUNTS',overlay=False).add_to(m) # 행정동별 총 자전거 반입량
# feature_group3 = folium.FeatureGroup(name='STATION_LEAS_COUNTS',overlay=False).add_to(m) # 행정동별 총 자전거 반출량
# feature_group4 = folium.FeatureGroup(name='RIDE_PER_POP',overlay=False).add_to(m) # 행정동별 한명당 자전거 이용량
fs = [feature_group1]#,feature_group1,feature_group2,feature_group3,feature_group4]
feature_list = ["STATION_RTN-LEAS_MEAN"]#,"STATION_RTN-LEAS_MEAN","STATION_RTN_COUNTS","STATION_LEAS_COUNTS","RIDE_PER_POP"]
for i in range(len(feature_list)):
choropleth1 = folium.Choropleth(
geo_data=gdf_boundary_town,
name='choropleth',
data=gdf_boundary_town,
columns=['DONG_NM', feature_list[i]],
key_on='feature.properties.DONG_NM',
fill_color='YlGn',
nan_fill_color="black",
fill_opacity=0.7,
line_opacity=0.2,
legend_name='Arrival (in Quintal)',
highlight=True,
line_color='black').geojson.add_to(fs[i])
colormap = linear.YlGn_09.scale(
gdf_boundary_town[feature_list[i]].min(),
gdf_boundary_town[feature_list[i]].max()).to_step(10)
colormap.caption = feature_list[i]
colormap.add_to(m)
folium.TileLayer('cartodbpositron',overlay=True,name="background").add_to(m)
folium.LayerControl('topright',collapsed=False).add_to(m)
m
결과적으로, 거주지역과 상업지역에 따른 반입, 반출의 특징이 상이한 것을 확인할 수 있었다.
2.1.1.2절에서 진행한 상업/거주지역에 따른 이용량을 비교해 본다.
| variable | means |
|---|---|
| TOTAL | 절대 수요(이용량) |
| pop_200_sum | 200m 내 주변 인구 대비 이용량 |
| RTN-LEAS | 상대 수요(이용량) |
gdf_station_cover["area_station_type"] = gdf_station_cover["HOUSE_COUNTS_200_mean"].apply(lambda x: 1 if x <= 1.0 else 0)
gdf_station_cover.tail()
sns.set(style="white", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=1,ncols=3)
fig.set_size_inches(12, 5)
sns.boxplot(x="area_station_type", y="TOTAL", data=gdf_station_cover, ax=axes[0])
sns.boxplot(x="area_station_type", y="pop_200_sum", data=gdf_station_cover, ax=axes[1])
sns.boxplot(x="area_station_type", y="RTN-LEAS", data=gdf_station_cover, ax=axes[2])
plt.tight_layout()
plt.show()
위의 boxplot을 가지고 알 수 있는 점은 다음과 같다.
정류장의 고도 분포를 확인하고, 상한 하한을 설정한다.
sns.distplot(station_df["altitude"])
bike_station_mean_altitude = np.mean(station_df["altitude"])
print("현재 정류장의 평균 고도 : ", bike_station_mean_altitude)
bike_station_std_altitude = np.std(station_df["altitude"])
print("현재 정류장의 표준편차 : ", bike_station_std_altitude)
굉장히 skew한 분포임을 알 수 있다.
10~20m에 많은 정류장이 분포하고 있음을 알 수 있다.
# 정류장 고도 기반 total 수요 보기
gdf_station_altitude = gdf_station[['altitude', 'TOTAL']]
gdf_station_altitude.describe().transpose()
# 고도에 따른 수요 plot
bound = 25
gdf_station_altitude['altitude_bounds'] = pd.cut(gdf_station['altitude'], bins=[8, bound, 44], labels=['low', 'high'])
gdf_station_altitude.head()
# 1분위수에 해당하는 25m를 기준으로 group을 나누고 boxplot을 그려보았다.
sns.boxplot(x="altitude_bounds",
y="TOTAL",
data=gdf_station_altitude)
plt.show()
print("LOW 평균 수요 : " , gdf_station_altitude[gdf_station_altitude['altitude_bounds'] == 'low']['TOTAL'].mean())
print("HIGH 평균 수요 : " , gdf_station_altitude[gdf_station_altitude['altitude_bounds'] == 'high']['TOTAL'].mean())
print("LOW 갯수 : " , gdf_station_altitude[gdf_station_altitude['altitude_bounds'] == 'low']['TOTAL'].count())
print("HIGH 갯수 : " , gdf_station_altitude[gdf_station_altitude['altitude_bounds'] == 'high']['TOTAL'].count())
1분위 수에 해당하는 25m를 기준으로 boxplot을 그려보았을 때, 두 집단의 평균차이가 유의미하게 남을 확인할 수 있다.
이를 기반으로 접근성이 용이한 기준을 약 25m로 선정할 것이다.
자전거 거치대 주변에 있는 대중교통과의 연관성을 파악한다.
2.3.2.1 근처 버스정류장 수 2.3.2.2 근처 버스정류장 승하차 평균 2.3.2.3 근처 지하철 수 2.3.2.4 근처 지하철 승하차 평균
import copy
import time
from haversine import haversine
def countCloseLandmarks(df1, df2, boundary, prefix):
global check
"""
<input>
df1 : 주변 정보를 셀 위치 정보(자전거 정류장)
df2 : 주변 위치 정보(버스정류장, 지하철역))
boundary :
prefix : 생성될 변수명 앞에 붙을 접두어
"""
# 1. df2의 lat, lon을 str로 변환한다.
# 2. 이걸 기반으로 lon_lat 변수를 만든다.
# 3. 거리를 계산한 np.array를 반환한다.
# 4. 그 np.array에서 조건을 걸고 넘는 것의 개수를 반환하여 ret_df에 추가한다.
ret_df = copy.deepcopy(df1)
count_df = copy.deepcopy(df2)
count_df["LON"] = count_df["LON"].apply(lambda x: str(x))
count_df["LAT"] = count_df["LAT"].apply(lambda x: str(x))
count_df["LAT_LON"] = count_df["LAT"] + ' ' + count_df["LON"]
count_loc_list = count_df["LAT_LON"].to_numpy()
start = time.time()
for i in range(ret_df.shape[0]):
if ret_df.shape[0] > 1000 and i % 1000 == 0:
print("진행중 {} / {} ".format(i , ret_df.shape[0]))
print("진행 시간", time.time() - start)
tar_loc = (ret_df.loc[i, "LAT"], ret_df.loc[i, "LON"])
f = lambda x: haversine(tar_loc, (float(x.split(' ')[0]), float(x.split(' ')[1])), unit = 'km')
result = np.fromiter((f(count_loc) for count_loc in count_loc_list), np.float, count=len(count_loc_list))
check = result
cnt_loc = sum(result <= boundary)
ret_df.loc[i, f'{prefix}_cnt_{int(boundary*1000)}m'] = cnt_loc
return ret_df
import copy
def calculateCloseLandmarkFeatures(df1, df2, boundary, prefix, feature, method=sum):
"""
<describe>
feature에 있는 nan를 지운 상태에서 승하차 인원 수를 더한다.
<input>
df1 : 주변 정보를 셀 위치 정보(자전거 정류장)
df2 : 주변 위치 정보(버스정류장, 지하철역))
boundary :
prefix : 생성될 변수명 앞에 붙을 접두어
feature : 연산을 수행할 feature (승하차)
method : sum, np.mean
"""
# 1. df2의 lat, lon을 str로 변환한다.
# 2. 이걸 기반으로 lon_lat 변수를 만든다.
# 3. 거리를 계산한 np.array를 반환한다.
# 4. 그 np.array에서 조건을 걸고 넘는 것의 데이터를 검색하여 연산을 진행한다.
ret_df = copy.deepcopy(df1)
count_df = copy.deepcopy(df2)
count_df = count_df[count_df[feature].notnull()]
count_df["LON"] = count_df["LON"].apply(lambda x: str(x))
count_df["LAT"] = count_df["LAT"].apply(lambda x: str(x))
count_df["LAT_LON"] = count_df["LAT"] + ' ' + count_df["LON"]
count_loc_list = count_df["LAT_LON"].to_numpy()
start = time.time()
for i in range(ret_df.shape[0]):
if ret_df.shape[0] > 1000 and i % 1000 == 0:
print("진행중 {} / {} ".format(i , ret_df.shape[0]))
print("진행 시간", time.time() - start)
tar_loc = (ret_df.loc[i, "LAT"], ret_df.loc[i, "LON"])
f = lambda x: haversine(tar_loc, (float(x.split(' ')[0]), float(x.split(' ')[1])), unit = 'km')
result = np.fromiter((f(count_loc) for count_loc in count_loc_list), np.float, count=len(count_loc_list))
cal_result = method(count_df.loc[(result <= boundary), feature])
ret_df.loc[i, f'{prefix}_{feature}_{method.__name__}_{int(boundary*1000)}m'] = cal_result
return ret_df
station_eda_df = copy.deepcopy(station_df)
station_eda_df = countCloseLandmarks(station_eda_df, busstop_df, 0.2, "busstop")
station_eda_df = countCloseLandmarks(station_eda_df, busstop_df, 0.3, "busstop")
station_eda_df = countCloseLandmarks(station_eda_df, busstop_df, 0.5, "busstop")
station_eda_df.head()
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=1,ncols=3)
fig.set_size_inches(12, 5)
x200 = station_eda_df["busstop_cnt_200m"]
x300 = station_eda_df["busstop_cnt_300m"]
x500 = station_eda_df["busstop_cnt_500m"]
sns.distplot(x200, color="m", ax=axes[0])
sns.distplot(x300, color="b", ax=axes[1])
sns.distplot(x500, color="r", ax=axes[2])
axes[0].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
axes[1].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
axes[2].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# plt.setp(axes, yticks=[])
plt.tight_layout()
plt.show()
200, 300, 500m 단위로 정류장 수를 시각화해보았을 때, 200, 300은 약간의 왜도를 가지는 분포가 생성되었으며, 500m를 보면 정규분포에 근사한 것을 확인할 수 있다.
200m기준으로 약 2~5개의 정류장이 분포한 위치가 사용량이 많다.
station_eda_df = copy.deepcopy(station_df)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, busstop_df, 0.2, "busstop", "GETON_CNT", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, busstop_df, 0.3, "busstop", "GETON_CNT", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, busstop_df, 0.5, "busstop", "GETON_CNT", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, busstop_df, 0.2, "busstop", "GETON_CNT", np.mean)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, busstop_df, 0.3, "busstop", "GETON_CNT", np.mean)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, busstop_df, 0.5, "busstop", "GETON_CNT", np.mean)
station_eda_df.head()
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=2,ncols=3)
fig.set_size_inches(12, 10)
x200 = station_eda_df["busstop_GETON_CNT_sum_200m"]
x300 = station_eda_df["busstop_GETON_CNT_sum_300m"]
x500 = station_eda_df["busstop_GETON_CNT_sum_500m"]
a = sns.distplot(x200, color="m", ax=axes[0,0])
b = sns.distplot(x300, color="b", ax=axes[0,1])
c = sns.distplot(x500, color="r", ax=axes[0,2])
axes[0,0].set_xticks([0, 200000, 400000, 600000, 800000])
axes[0,1].set_xticks([0, 200000, 400000, 600000, 800000])
axes[0,2].set_xticks([0, 200000, 400000, 600000, 800000])
x200 = station_eda_df["busstop_GETON_CNT_mean_200m"]
x300 = station_eda_df["busstop_GETON_CNT_mean_300m"]
x500 = station_eda_df["busstop_GETON_CNT_mean_500m"]
a = sns.distplot(x200, color="m", ax=axes[1,0])
b = sns.distplot(x300, color="b", ax=axes[1,1])
c = sns.distplot(x500, color="r", ax=axes[1,2])
# a.set_xticklabels(rotation=30)
# axes[0].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[1].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[2].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# plt.setp(axes, yticks=[])
plt.tight_layout()
plt.show()
굉장히 왜도가 높은 분포가 형성되었다.
이는 대부분의 자전거 정류장이 사용량이 많지 않은 정류장을 대거 포함하고 있다는 사실을 대변한다.
# 지하철 관련
station_eda_df = copy.deepcopy(station_df)
station_eda_df = countCloseLandmarks(station_eda_df, subway_space_df, 0.5, "subway")
station_eda_df = countCloseLandmarks(station_eda_df, subway_space_df, 1.0, "subway")
station_eda_df = countCloseLandmarks(station_eda_df, subway_space_df, 1.5, "subway")
station_eda_df.head()
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=1,ncols=3)
fig.set_size_inches(12, 5)
x200 = station_eda_df["subway_cnt_500m"]
x300 = station_eda_df["subway_cnt_1000m"]
x500 = station_eda_df["subway_cnt_1500m"]
sns.distplot(x200, color="m", ax=axes[0])
sns.distplot(x300, color="b", kde=False, ax=axes[1])
sns.distplot(x500, color="r", ax=axes[2])
# axes[0].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[1].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[2].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# plt.setp(axes, yticks=[])
plt.tight_layout()
plt.show()
버스 정류장과 다르게, 자전거 정류장 주변에 지하철 역은 많이 분포하지 않는 것으로 보인다.
station_eda_df = copy.deepcopy(station_df)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 0.5, "subway", "IN", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.0, "subway", "IN", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.5, "subway", "IN", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 0.5, "subway", "IN", np.mean)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.0, "subway", "IN", np.mean)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.5, "subway", "IN", np.mean)
station_eda_df.head()
# 지하철역 승차인원 시각화
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=2,ncols=3)
fig.set_size_inches(12, 10)
x200 = station_eda_df["subway_IN_sum_500m"]
x300 = station_eda_df["subway_IN_sum_1000m"]
x500 = station_eda_df["subway_IN_sum_1500m"]
sns.distplot(x200, color="m", ax=axes[0,0])
sns.distplot(x300, color="b", ax=axes[0,1])
sns.distplot(x500, color="r", ax=axes[0,2])
x200 = station_eda_df["subway_IN_mean_500m"]
x300 = station_eda_df["subway_IN_mean_1000m"]
x500 = station_eda_df["subway_IN_mean_1500m"]
sns.distplot(x200, color="m", ax=axes[1,0])
sns.distplot(x300, color="b", ax=axes[1,1])
sns.distplot(x500, color="r", ax=axes[1,2])
plt.tight_layout()
plt.show()
반경을 넓힐 수록 정규분포에 근사함을 확인할 수 있었다.
다만 반경을 작게 했을 때, 승차 인원이 작은 부분에 분포가 몰렸다.
이는 승차인원이 많은 큰 지하철 역은 보다 멀리 자전거 정류장이 위치해 있다는 사실을 반증한다.
station_eda_df = copy.deepcopy(station_df)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 0.5, "subway", "OUT", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.0, "subway", "OUT", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.5, "subway", "OUT", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 0.5, "subway", "OUT", np.mean)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.0, "subway", "OUT", np.mean)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.5, "subway", "OUT", np.mean)
station_eda_df.head()
# 지하철역 하차인원 시각화
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=2,ncols=3)
fig.set_size_inches(12, 10)
x200 = station_eda_df["subway_OUT_sum_500m"]
x300 = station_eda_df["subway_OUT_sum_1000m"]
x500 = station_eda_df["subway_OUT_sum_1500m"]
sns.distplot(x200, color="m", ax=axes[0,0])
sns.distplot(x300, color="b", ax=axes[0,1])
sns.distplot(x500, color="r", ax=axes[0,2])
x200 = station_eda_df["subway_OUT_mean_500m"]
x300 = station_eda_df["subway_OUT_mean_1000m"]
x500 = station_eda_df["subway_OUT_mean_1500m"]
sns.distplot(x200, color="m", ax=axes[1,0])
sns.distplot(x300, color="b", ax=axes[1,1])
sns.distplot(x500, color="r", ax=axes[1,2])
plt.tight_layout()
plt.show()
하차 인원도 승차 인원과 같은 양상을 보인다.
station_eda_df = copy.deepcopy(station_df)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 0.5, "subway", "NET_TOTAL", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.0, "subway", "NET_TOTAL", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.5, "subway", "NET_TOTAL", sum)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 0.5, "subway", "NET_TOTAL", np.mean)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.0, "subway", "NET_TOTAL", np.mean)
station_eda_df = calculateCloseLandmarkFeatures(station_eda_df, subway_space_df, 1.5, "subway", "NET_TOTAL", np.mean)
station_eda_df.head()
# 상대 수요
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=2,ncols=3)
fig.set_size_inches(12, 10)
x200 = station_eda_df["subway_NET_TOTAL_sum_500m"]
x300 = station_eda_df["subway_NET_TOTAL_sum_1000m"]
x500 = station_eda_df["subway_NET_TOTAL_sum_1500m"]
sns.distplot(x200, color="m", ax=axes[0,0])
sns.distplot(x300, color="b", ax=axes[0,1])
sns.distplot(x500, color="r", ax=axes[0,2])
x200 = station_eda_df["subway_NET_TOTAL_mean_500m"]
x300 = station_eda_df["subway_NET_TOTAL_mean_1000m"]
x500 = station_eda_df["subway_NET_TOTAL_mean_1500m"]
sns.distplot(x200, color="m", ax=axes[1,0])
sns.distplot(x300, color="b", ax=axes[1,1])
sns.distplot(x500, color="r", ax=axes[1,2])
plt.tight_layout()
plt.show()
총 유입량에 대한 분포를 그려보았을 때, 대부분 0 근처에 높은 값이 산출되었다.
이는 대부분의 지하철에서 순 사용량이 비슷하다는 것을 말한다.
하지만 양 극단에서는 유립량이 +, -가 발생하는 것으로 보아, 이러한 지하철 역의 특성을 반영해주어야 하는 필요성을 알 수 있다.
지하철 역에 대한 분포를 보았을 때, 버스 정류장에 비해 접근성이 비교적 떨어진다는 결론을 얻을 수 있었다.
하지만, 이러한 역의 존재는 자전거 정류장에 분명한 영향력을 줄 수 있으므로
유입량이 많은 역, 적은 역과 같은 특성을 모델 훈련에 있어 사용해야 한다는 결론을 얻었다.
candidate = ['초등학교', '중학교', '고등학교', '특수학교', '각종학교', '대학']
str_expr = f"TYPE in {candidate}"
school_gdf = city_plan_public_gdf.query(str_expr) # 조건 부합 데이터 추출
school_gdf.reset_index(drop=True, inplace=True)
print(school_gdf.shape)
school_gdf.head()
school_gdf["CENTROID"] = school_gdf["geometry"].apply(lambda x: x.centroid)
school_gdf["LON"] = school_gdf["CENTROID"].apply(lambda x: x.x)
school_gdf["LAT"] = school_gdf["CENTROID"].apply(lambda x: x.y)
print(school_gdf.shape)
school_gdf.head()
# 교육시설 관련
station_eda_df = copy.deepcopy(station_df)
station_eda_df = countCloseLandmarks(station_eda_df, school_gdf, 0.5, "school")
station_eda_df = countCloseLandmarks(station_eda_df, school_gdf, 1.0, "school")
station_eda_df = countCloseLandmarks(station_eda_df, school_gdf, 1.5, "school")
station_eda_df.head()
# 거치대 주변 교육시설 시각화
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=1,ncols=3)
fig.set_size_inches(12, 5)
x200 = station_eda_df["school_cnt_500m"]
x300 = station_eda_df["school_cnt_1000m"]
x500 = station_eda_df["school_cnt_1500m"]
sns.distplot(x200, color="m", kde=False, ax=axes[0])
sns.distplot(x300, color="b", ax=axes[1])
sns.distplot(x500, color="r", ax=axes[2])
# axes[0].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[1].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[2].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# plt.setp(axes, yticks=[])
plt.tight_layout()
plt.show()
학교 같은 경우에는 정류장 주변 개수가 200m 기준 적어도 1개정도 이상을 가지고 있다.
이는 매우 유의미한 결과이다.
concerthall_space_df_temp = concerthall_space_df[["NAME", "LAT", "LON"]]
museum_space_df_temp = museum_space_df[["NAME", "LAT", "LON"]]
culture_facility_df = concerthall_space_df_temp.append(museum_space_df_temp)
culture_facility_df.reset_index(drop=True, inplace=True)
culture_facility_df.head()
# 문화시설 관련
station_eda_df = copy.deepcopy(station_df)
station_eda_df = countCloseLandmarks(station_eda_df, culture_facility_df, 0.5, "cul_fac")
station_eda_df = countCloseLandmarks(station_eda_df, culture_facility_df, 1.0, "cul_fac")
station_eda_df = countCloseLandmarks(station_eda_df, culture_facility_df, 1.5, "cul_fac")
station_eda_df.head()
# 거치대 주변 문화시설 시각화
import seaborn as sns, numpy as np
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=1,ncols=3)
fig.set_size_inches(12, 5)
x200 = station_eda_df["cul_fac_cnt_500m"]
x300 = station_eda_df["cul_fac_cnt_1000m"]
x500 = station_eda_df["cul_fac_cnt_1500m"]
sns.distplot(x200, color="m", kde=False, ax=axes[0])
sns.distplot(x300, color="b", ax=axes[1])
sns.distplot(x500, color="r", ax=axes[2])
# axes[0].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[1].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[2].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# plt.setp(axes, yticks=[])
plt.tight_layout()
plt.show()
문화시설의 개수를 500, 1000, 1500m 기준으로 보았을 때, 대부분이 0개를 가지고 있는 것을 알 수 있다.
# 체육시설 관련
station_eda_df = copy.deepcopy(station_df)
station_eda_df = countCloseLandmarks(station_eda_df, sport_complex_info_df, 0.5, "sports")
station_eda_df = countCloseLandmarks(station_eda_df, sport_complex_info_df, 1.0, "sports")
station_eda_df = countCloseLandmarks(station_eda_df, sport_complex_info_df, 1.5, "sports")
station_eda_df.head()
# 주변 체육시설 시각화
sns.set(style="whitegrid", palette="muted", color_codes=True)
fig, axes = plt.subplots(nrows=1,ncols=3)
fig.set_size_inches(12, 5)
x200 = station_eda_df["sports_cnt_500m"]
x300 = station_eda_df["sports_cnt_1000m"]
x500 = station_eda_df["sports_cnt_1500m"]
sns.distplot(x200, color="m", kde=False, ax=axes[0])
sns.distplot(x300, color="b", ax=axes[1])
sns.distplot(x500, color="r", ax=axes[2])
# axes[0].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[1].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# axes[2].set_yticks([0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3])
# plt.setp(axes, yticks=[])
plt.tight_layout()
plt.show()
체육 시설 역시 500m 기준 70% 이상이 0이다.
결과적으로 대부분의 문화시설과 체육시설이 이를 기준으로 배치되지 않았음을 확인할 수 있었다.
오히려 학교와 같은 교육시설을 기반으로 많은 수의 정류장이 배치되었음을 확인했다.
도로를 구분하는 다양한 구분 코드를 기반으로 시각화하면서, 어떤 특징이 있는지 확인한다.
| code | meaning |
|---|---|
| 1 | 행자부도로 |
| 2 | 광역도로 |
| 3 | 시군구도로 |
| code | meaning |
|---|---|
| 0 | 주도로 |
| 1 | 1차 종속도로 |
| 2 | 2차 종속도로 |
| code | meaning |
|---|---|
| 1 | 고속도로 |
| 2 | 대로 |
| 3 | 로 |
| 4 | 길 |
# folium marker를 html로 변환하는 코드
table = """
<!DOCTYPE html>
<html>
<head>
<style>
table {{
width:100%;
}}
table, th, td {{
border: 1px solid black;
border-collapse: collapse;
}}
th, td {{
padding: 5px;
text-align: left;
}}
table#t01 tr:nth-child(odd) {{
background-color: #eee;
}}
table#t01 tr:nth-child(even) {{
background-color:#fff;
}}
</style>
</head>
<body>
<table id="t01">
<tr>
<td>STATION ID</td>
<td>{}</td>
</tr>
<tr>
<td>Name</td>
<td>{}</td>
</tr>
<tr>
<td>Holder Capacity</td>
<td>{}</td>
</tr>
<tr>
<td>Altitude</td>
<td>{}</td>
</tr>
</table>
</body>
</html>
""".format
####################### 해당 환경에서 제대로 동작하지 않음 ##############################
import random as rd
import folium
from folium import IFrame
from folium.plugins import MarkerCluster, FeatureGroupSubGroup
point_num = 100
# 다 제작했으나 해당 환경에 구동하지 않음
# 광역도로 구분 코드(WDR_RD_CD) road_code1_df
# 도로구간종속구분코드(RDS_DPN_SE) road_code2_df
# 도로위계구분코드(ROA_CLS_SE) road_code3_df
m = folium.Map(location=center, tiles="OpenStreetMap", zoom_start=12)
fg = folium.FeatureGroup(name='ALL')
m.add_child(fg)
g0 = FeatureGroupSubGroup(fg, 'Bike station[star]')
m.add_child(g0)
g1 = FeatureGroupSubGroup(fg, 'Gwangyeok 1 - Hangjaboo[skyblue]')
m.add_child(g1)
g2 = FeatureGroupSubGroup(fg, 'Gwangyeok 2 - Gwangyeok[brue]')
m.add_child(g2)
g3 = FeatureGroupSubGroup(fg, 'Gwangyeok 3 - SiGunGu[red]')
m.add_child(g3)
g4 = FeatureGroupSubGroup(fg, 'Subjected Road 0 - Main[green]')
m.add_child(g4)
g5 = FeatureGroupSubGroup(fg, 'Subjected Road 1 - 1st Subjected[springgreen]')
m.add_child(g5)
g6 = FeatureGroupSubGroup(fg, 'Subjected Road 1 - 2nd Subjected[violet]')
m.add_child(g6)
g7 = FeatureGroupSubGroup(fg, 'Road Level 1 - Highway[grey]')
m.add_child(g7)
g8 = FeatureGroupSubGroup(fg, 'Road Level 2 - Dae Ro[pink]')
m.add_child(g8)
g9 = FeatureGroupSubGroup(fg, 'Road Level 3 - Ro[lime]')
m.add_child(g9)
g10 = FeatureGroupSubGroup(fg, 'Road Level 4 - Gil[yellow]')
m.add_child(g10)
width, height = 310,200
# 버스정류장 위치
popups, locations = [], []
# marker_cluster = MarkerCluster().add_to(g0)
# for idx, row in station_df.iterrows():
# folium.Marker(
# location=[row['geometry'].y, row['geometry'].x],
# popup=IFrame(table(row['STATION_ID'], row['STATION_NAME'], row['HOLDER_CAPACITY'], row['altitude']), width=width, height=height),
# icon=folium.Icon(color='red',icon='star'),
# ).add_to(marker_cluster)
for idx, row in station_df.iterrows():
locations.append([row['geometry'].y, row['geometry'].x])
num = row['STATION_ID']
name = row['STATION_NM']
holder = row["HOLDER_CAPACITY"]
altitude = row["altitude"]
iframe = IFrame(table(num, name, holder, altitude), width=width, height=height)
popups.append(iframe)
g0.add_children(MarkerCluster(locations=locations, popups=popups))
m.add_children(g0)
# 광역도로
# folium.GeoJson(
# address_road_gdf[address_road_gdf["WDR_RD_CD"] == "1"],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "skyblue",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g1)
# folium.GeoJson(
# address_road_gdf[address_road_gdf["WDR_RD_CD"] == "2"],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "brue",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g2)
folium.GeoJson(
address_road_gdf.loc[rd.sample(list(address_road_gdf[address_road_gdf["WDR_RD_CD"] == "3"].index), point_num)],
style_function=lambda feature: {
'fillOpacity': 0.1,
'weight': 3,
"color": "black",
'lineColor': '#E2E2E2',
"opacity": 0.3,
}
).add_to(g3)
# # 종속구간
# folium.GeoJson(
# address_road_gdf.loc[rd.sample(list(address_road_gdf[address_road_gdf["RDS_DPN_SE"] == "0"].index), point_num)],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "green",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g4)
# folium.GeoJson(
# address_road_gdf[address_road_gdf["RDS_DPN_SE"] == "1"],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "springgreen",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g5)
# folium.GeoJson(
# address_road_gdf[address_road_gdf["RDS_DPN_SE"] == "2"],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "violet",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g6)
# # 위계코드
# folium.GeoJson(
# address_road_gdf[address_road_gdf["ROA_CLS_SE"] == "1"],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "grey",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g7)
# folium.GeoJson(
# address_road_gdf[address_road_gdf["ROA_CLS_SE"] == "2"],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "pink",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g8)
# folium.GeoJson(
# address_road_gdf[address_road_gdf["ROA_CLS_SE"] == "3"],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "lime",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g9)
# folium.GeoJson(
# address_road_gdf[address_road_gdf["ROA_CLS_SE"] == "4"],
# style_function=lambda feature: {
# 'fillOpacity': 0.1,
# 'weight': 3,
# "color": "yellow",
# 'lineColor': '#E2E2E2',
# "opacity": 0.3,
# }
# ).add_to(g10)
folium.GeoJson(
administration_boundary_town_gdf,
name='Administration Boundary Layer',
style_function=lambda feature: {
'fillOpacity': 0.2,
'weight': 2,
"color": "#1C1C1C",
'lineColor': '#E2E2E2',
"opacity": 0.4,
}
).add_to(m)
folium.LayerControl(collapsed=False).add_to(m)
m
folium 렌더링의 오류가 있어 해당 데이터를 100개만 시각화해보았다.
모든 도로 구분 코드에 대해 분류하여 살펴본 결과, 시군구 도로(광역 코드 3번)는 자전거를 타기 용이한 모든 도로에 대해 표현됨을 확인했다.
공공 자전거 공유 문제에서는 필연적으로 수요-공급 불균형 문제가 발생한다. 이러한 점을 추후에 해결하기 위해, 어떠한 특징을 갖는지 분석한다.
demand_month_df = createDemandDf(route_df, unit='month')
demand_month_df
demand_month = demand_month_df[['MONTH', 'LEAS_COUNTS', 'RTN_COUNTS']].groupby('MONTH').mean()
demand_month
demand_month["month"] = demand_month.index
demand_month
import seaborn as sns
f = plt.figure(figsize=(10, 6))
gs = f.add_gridspec(1, 2)
with sns.axes_style("whitegrid"):
ax = f.add_subplot(gs[0, 0])
sns.barplot(x="month", y="LEAS_COUNTS", data=demand_month, palette="Greens", ax=ax)
with sns.axes_style("whitegrid"):
ax = f.add_subplot(gs[0, 1])
sns.barplot(x="month", y="RTN_COUNTS", data=demand_month, palette="Blues", ax=ax)
demand_day_df = createDemandDf(route_df, unit='day')
demand_day_df
demand_day = demand_day_df[['DAY', 'LEAS_COUNTS', 'RTN_COUNTS']].groupby('DAY').mean()
demand_day
demand_day["day"] = demand_day.index
demand_day
import seaborn as sns
f = plt.figure(figsize=(10, 6))
gs = f.add_gridspec(1, 2)
with sns.axes_style("whitegrid"):
ax = f.add_subplot(gs[0, 0])
sns.barplot(x="day", y="LEAS_COUNTS", data=demand_day, palette="Greens", ax=ax)
with sns.axes_style("whitegrid"):
ax = f.add_subplot(gs[0, 1])
sns.barplot(x="day", y="RTN_COUNTS", data=demand_day, palette="Blues", ax=ax)
demand_date_df = createDemandDf(route_df)
demand_date_df
demand_date_df["DATE"] = demand_date_df["DATE"].apply(lambda x: int(x.day))
demand_date_df
demand_date = demand_date_df[['DATE', 'LEAS_COUNTS', 'RTN_COUNTS']].groupby('DATE').mean()
demand_date.head()
demand_date["date"] = demand_date.index
demand_date.head()
import seaborn as sns
f = plt.figure(figsize=(10, 6))
gs = f.add_gridspec(1, 2)
with sns.axes_style("whitegrid"):
ax = f.add_subplot(gs[0, 0])
sns.barplot(x="date", y="LEAS_COUNTS", data=demand_date, palette="Greens", ax=ax)
with sns.axes_style("whitegrid"):
ax = f.add_subplot(gs[0, 1])
sns.barplot(x="date", y="RTN_COUNTS", data=demand_date, palette="Blues", ax=ax)
route_df['LEA_HOUR'] = route_df['LEAS_DATE'].apply(lambda x : x.hour)
route_df['RTN_HOUR'] = route_df['RTN_DATE'].apply(lambda x : x.hour)
lea_hour = route_df[['LEA_HOUR', 'LEAS_NO']].groupby('LEA_HOUR').count()
rtn_hour = route_df[['RTN_HOUR', 'LEAS_NO']].groupby('RTN_HOUR').count()
lea_hour["HOUR"] = lea_hour.index
rtn_hour["HOUR"] = rtn_hour.index
lea_hour.columns = ['LEAS_COUNTS', 'HOUR']
rtn_hour.columns = ['RTN_COUNTS', 'HOUR']
import seaborn as sns
f = plt.figure(figsize=(13, 6))
gs = f.add_gridspec(1, 2)
with sns.axes_style("whitegrid"):
ax = f.add_subplot(gs[0, 0])
sns.barplot(x="HOUR", y="LEAS_COUNTS", data=lea_hour, palette="Greens", ax=ax)
with sns.axes_style("whitegrid"):
ax = f.add_subplot(gs[0, 1])
sns.barplot(x="HOUR", y="RTN_COUNTS", data=rtn_hour, palette="Blues", ax=ax)
위의 네가지 결과를 보고 알 수 있는 점은 다음과 같다.
시간대 별 차이는 조금더 확대해서 본다면 출근, 퇴근 시간과 연관지어 생각할 수 있다. 또한 이런 출퇴근은, 거주지역, 상업지역에 따라 변화할 수 있다.
각각의 거치대 이용량의 일별 평균(상대수요)을 기반으로 거치대 별 반입/반출에 대해 영향력을 파악한다.
sns.scatterplot(x="STATION_ID", y="RTN-LEAS_MEAN", data=gdf_station[:-1])
대체적으로 거치대에서 반출이 많은 것을 파악할 수 있다. 전체 자전거양이 유지된다는 가정을 도입한다면, 이러한 분포에는 문제가 있다는 것을 파악할 수 있다.
print("상대 수요 평균 : ", gdf_station['RTN-LEAS_MEAN'][:-1].mean())
앞서 분석을 진행하며, 대부분의 거치대에서 반출이 많다는 사실을 파악했다. 거치대 중 특이한 것이 있는지 분석하고, 이를 어떻게 활용할지 파악한다.
앞서 상대 수요 평균이 −2.3 이 나오는 것은 모든 거치대에서 반출이 평균적으로 -2.3개 정도 많다고 생각할 수 있다. 하지만 자전거의 총량이 변화하지 않는다는 가정을 도입한다면 이는 현실적으로 가능하지 않다. 따라서 0번 station의 의미를 파악해 본다.
print("LEAS_STATION 이 0번 스테이션인 경우")
route_df[route_df['LEAS_STATION'] == 0].shape
print("RTN_STATION 이 0번 스테이션인 경우")
route_df[route_df['RTN_STATION'] == 0].shape
0번 스테이션은 반출은 없고 반납만 존재하는 station이다.
station_num = 155
years = 365
print("0번 스테이션의 평균 : ", len(route_df[route_df['RTN_STATION'] == 0]) / station_num / (years * 3))
station_num = 155
years = 365
print("0번 스테이션의 평균 : ", len(route_df[route_df['RTN_STATION'] == 0]) / station_num / (years * 3))
이로 부터 알 수 있는 사실은 다음과 같다.
이러한 점을 근거로 우리는 이 0번 거치대는 중앙에서 거치대의 Rebalancing을 해소하기 위한 운반이라고 가정하였다. 이러한 점을 가정한다면, 0번 거치대 이외의 모든 거치대의 상대수요는 평균적으로 2.3 더해져야 할 것이다.
bias = 2.3
import copy
gdf_station_temp_visualize = copy.deepcopy(gdf_station)
gdf_station_temp_visualize["RTN-LEAS_MEAN"] = gdf_station_temp_visualize["RTN-LEAS_MEAN"].apply(lambda x: x+bias)
sns.scatterplot(x="STATION_ID", y="RTN-LEAS_MEAN", data=gdf_station_temp_visualize[:-1])
bias가 도입된 산점도는 위와 같다.
추후 이러한 가정을 도입하여 모델링을 진행할 것이다.
거치대간 상호작용을 네트워크 분석을 통해 파악한다.
# 자전거 정류장 간 network counting하는 함수
def createNetworkDf(route_df, unit='all', remove_stations=[0, 998, 999]):
'''
input // route_df = 01.운영이력.csv, unit = [all, month, day] // day = 0 => 월요일
output // columns = LEAS_STATION, RTN_STATION, COUNTS // form : pd.DataFrame
time // unit = 'all' ==> 0.2sec , unit = 'month' ==> 32sec , unit = 'day' ==> 30sec
'''
possible_unit = ['all', 'month', 'day']
if unit not in possible_unit:
print("unit Error!")
return
stations = sorted(route_df['LEAS_STATION'].unique())
for remove_station in remove_stations:
if remove_station in stations:
stations.remove(remove_station)
stations_num = len(stations)
if unit == 'all':
numpy_route_network = np.zeros((stations_num * stations_num, 2))
for i in range(stations_num):
for j in range(stations_num):
numpy_route_network[i*stations_num + j] = [stations[i], stations[j]]
route_network_df = pd.DataFrame(numpy_route_network, columns=['LEAS_STATION', 'RTN_STATION'])
temp_route_network_df = route_df[['LEAS_STATION', 'RTN_STATION', 'LEAS_NO']].groupby(['LEAS_STATION', 'RTN_STATION']).count()
temp_route_network_df.reset_index(inplace=True)
route_network_df = pd.merge(route_network_df, temp_route_network_df, how='left', on = ['LEAS_STATION', 'RTN_STATION'])
route_network_df.fillna(0, inplace=True)
return route_network_df
elif unit == 'month':
if 'MONTH' not in route_df.columns:
route_df['MONTH'] = route_df['LEAS_DATE'].apply(lambda x : datetime.strptime(x, "%Y-%m-%d %H:%M:%S").month)
numpy_route_network = np.zeros((12, stations_num * stations_num, 2))
for month in range(12):
for i in range(stations_num):
for j in range(stations_num):
numpy_route_network[month, i*stations_num + j] = [stations[i], stations[j]]
route_network_df = pd.DataFrame(columns=['LEAS_STATION', 'RTN_STATION', 'MONTH'])
for month in range(12):
month_network_df = pd.DataFrame(numpy_route_network[month], columns=['LEAS_STATION', 'RTN_STATION'])
month_network_df['MONTH'] = month + 1
route_network_df = route_network_df.append(month_network_df)
temp_route_network_df = route_df[['LEAS_STATION', 'RTN_STATION', 'MONTH', 'LEAS_NO']].groupby(['LEAS_STATION', 'RTN_STATION', 'MONTH']).count()
temp_route_network_df.reset_index(inplace=True)
route_network_df = pd.merge(route_network_df, temp_route_network_df, how='left', on = ['LEAS_STATION', 'RTN_STATION', 'MONTH'])
route_network_df.fillna(0, inplace=True)
return route_network_df
elif unit == 'day':
if 'DAY' not in route_df.columns:
route_df['DAY'] = route_df['LEAS_DATE'].apply(lambda x : datetime.strptime(x, "%Y-%m-%d %H:%M:%S").weekday())
numpy_route_network = np.zeros((7, stations_num * stations_num, 2))
for day in range(7):
for i in range(stations_num):
for j in range(stations_num):
numpy_route_network[day, i*stations_num + j] = [stations[i], stations[j]]
route_network_df = pd.DataFrame(columns=['LEAS_STATION', 'RTN_STATION', 'DAY'])
for day in range(7):
day_network_df = pd.DataFrame(numpy_route_network[day], columns=['LEAS_STATION', 'RTN_STATION'])
day_network_df['DAY'] = day
route_network_df = route_network_df.append(day_network_df)
temp_route_network_df = route_df[['LEAS_STATION', 'RTN_STATION', 'DAY', 'LEAS_NO']].groupby(['LEAS_STATION', 'RTN_STATION', 'DAY']).count()
temp_route_network_df.reset_index(inplace=True)
route_network_df = pd.merge(route_network_df, temp_route_network_df, how='left', on = ['LEAS_STATION', 'RTN_STATION', 'DAY'])
route_network_df.fillna(0, inplace=True)
return route_network_df
# 네트워크 분석에 필요한 데이터 프레임을 제작한다.
network_df = createNetworkDf(route_df, unit='all')
network_df['LEAS_STATION'] = network_df['LEAS_STATION'].apply(lambda x : int(x))
network_df['RTN_STATION'] = network_df['RTN_STATION'].apply(lambda x : int(x))
network_df = pd.merge(network_df, station_df[['STATION_ID', 'LAT', 'LON']], how='left', left_on=['LEAS_STATION'], right_on=['STATION_ID'])
network_df = pd.merge(network_df, station_df[['STATION_ID', 'LAT', 'LON']], how='left', left_on=['RTN_STATION'], right_on=['STATION_ID'])
network_df.columns = ['LEAS_STATION','RTN_STATION','COUNTS','STATION_ID_x','LEAS_LAT','LEAS_LONG','STATION_ID_y','RTN_LAT','RTN_LONG']
network_df = network_df[['LEAS_STATION','LEAS_LAT','LEAS_LONG','RTN_STATION','RTN_LAT','RTN_LONG','COUNTS']]
network_df['COUNTS'] = network_df['COUNTS'].apply(lambda x : int(x))
#network_df = network_df.loc[:1000]
network_df
# 간선과 가중치를 정의한다.
import networkx as nx
edges = network_df[['LEAS_STATION', 'RTN_STATION']].values.tolist()
weights = [int(l) for l in network_df['COUNTS'].values.tolist()]
G = nx.Graph(directed=True)
G.add_edges_from(edges)
for idx, node in enumerate(G.edges(data=True)):
G.edges[(node[0], node[1])]['weight'] = weights[idx] / network_df['COUNTS'].max()
import community as community_louvain
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import networkx as nx
# load the karate club graph
#G = nx.karate_club_graph()
# plt.figsize(20,20)
# 가장 좋은 클러스터(partition)으로 나눈다.
partition = community_louvain.best_partition(G)
# draw the graph
pos = nx.spring_layout(G)
# color the nodes according to their partition
cmap = cm.get_cmap('viridis', max(partition.values()) + 1)
nx.draw_networkx_nodes(G, pos, partition.keys(), node_size=40,
cmap=cmap, node_color=list(partition.values()))
nx.draw_networkx_edges(G, pos, alpha=0.1)
plt.show()
network_df.head()
단지 시각화만 해서는 어떠한 특징을 반영해 주기 어렵다.
우리는 이 네트워크 결과를 기반으로 가장 특성을 잘 나누는 6개의 군집으로 나누었다.
이를 추가적인 변수로 사용할 것이다.
# 네트워크 분석을 통해 도출된 최적 군집 6개에 대한 데이터 프레임을 제작한다.
cluster = pd.DataFrame([list(partition.keys()), list(partition.values())]).T
cluster.columns = ['STATION_ID', 'CLUSTER']
cluster
# 거치대 데이터에 해당 결과를 병합한다.
station_df = pd.merge(station_df, cluster, how='left', on='STATION_ID')
EDA 결과를 기반으로 예상되는 후보군을 추출한다.

사람들이 지나가는 도로를 시군구 도로로 제한하였다. 이는 EDA의 결과로 부터, 시군구 도로가 고양시의 주거지역, 상업지역 모두를 통과하는 Code라는 것을 확인했기 때문이다. 다른 조건을 기반으로 수행하지 않은 이유는, 시군구 도로가 고양시의 대부분의 도로의 특징을 모두 대변했기 때문이다. 가장 넓은 Coverage를 갖는 도로로부터 후보군을 추출함으로써 후보군의 일반화를 도모하고자 했다.
도로 주변 후보군을 고르는 간격은 300m를 기준으로 추출했다. 이 이유는 천천히 걸을 때 사람의 평균 속도를 시간당 3.5km라 가정했을 때, 5분 동안 갈 수 있는 거리이기 때문이다.
도로의 경계를 뽑기 위해서 도로의 폭을 고려했다. 제공된 데이터를 사용하고자 하였으나, 데이터의 문제가 있다는 것을 확인하고 기본적으로 1m를 기준으로 확장하여 판단하였다. 후에 도로 근처에 이를 옮김으로써 해결한다.
# WDR_RD_CD == 3 일때 모든 정류장 분포 (시군구 도로)
sigungu_road_gdf = address_road_gdf[address_road_gdf["WDR_RD_CD"] == "3"]
sigungu_road_gdf = sigungu_road_gdf.reset_index(drop=True)
sigungu_road_gdf.head()
# 도로 geometry 변환
sigungu_road_gdf = sigungu_road_gdf.set_crs(default_crs)
print("기본 crs :", sigungu_road_gdf.crs)
sigungu_road_meter_gdf = sigungu_road_gdf.to_crs(meter_crs)
print("변환된 crs : ", sigungu_road_meter_gdf.crs)
# 도로 -> Polygon화 -> boundary에서 300m 기반으로 점을 추출하는 함수
from shapely.geometry import Point
import math
def calcDis(pt1, pt2):
return Point(pt1).distance(Point(pt2))
#return np.sqrt((pt1[0]-pt2[0])**2 + (pt1[1]-pt2[1])**2)
def findInternal(pt1, pt2 ,dis):
if calcDis(pt1, pt2) < dis:
return (False, calcDis(pt1, pt2) - (dis - 300))
y = pt2[1] - pt1[1]
x = pt2[0] - pt1[0]
theta = math.atan2(y, x)
ny = pt1[1] + dis * math.sin(theta)
nx = pt1[0] + dis * math.cos(theta)
return (True, nx, ny)
def findCandidate(coords, dis=300):
ret_coords = []
#coords.append(coords[0])
coords_num = len(coords)
depot = dis
for i in range(coords_num - 2):
now = dis - depot
maxDis = calcDis(coords[i], coords[i + 1])
while True:
ret_val = findInternal(coords[i], coords[i + 1], now)
if ret_val[0] == False:
depot = ret_val[1]
#print(i, depot)
break
ret_coords.append((ret_val[1], ret_val[2]))
now += dis
depot = 0
# print(i)
return ret_coords
# 도로 -> 점 추출
buffer_size = 1
basic_type = type((sigungu_road_meter_gdf.loc[0, "geometry"]).buffer(buffer_size).boundary)
#candidate_gdf = gpd.GeoDataFrame(columns=["geometry"], crs=meter_crs)
num = len(sigungu_road_meter_gdf)
candidate_gdf = np.zeros((1, 2))
for idx, row in sigungu_road_meter_gdf.iterrows():
if idx % 1000 == 0:
print("진행률 : {} / {} ".format(idx, num ) )
road_data_meter = row["geometry"]
roadmap = road_data_meter.buffer(buffer_size).boundary
if type(roadmap) == basic_type:
coords = list(roadmap.coords)
temp = findCandidate(coords)
coords_list = np.array(temp)
candidate_gdf = np.concatenate((candidate_gdf, coords_list), axis=0)
else:
coords_num = len(roadmap)
for coords in roadmap:
coords = list(coords.coords)
temp = findCandidate(coords)
coords_list = np.array(temp)
candidate_gdf = np.concatenate((candidate_gdf, coords_list), axis=0)
candidate_gdf = list(map(Point, list(candidate_gdf[1:])))
candidate_gdf = {'geometry': candidate_gdf}
candidate_gdf = gpd.GeoDataFrame(candidate_gdf, crs=meter_crs)
candidate_gdf.head()
# 결과 확인
print(candidate_gdf.shape)
candidate_gdf.head()
# 후보군 plot
candidate_gdf.plot(markersize=1)
# 위경도 좌표계로 변환하기
candidate_gdf = candidate_gdf.to_crs(default_crs)
candidate_gdf
# 후보군 folium plot
# 노트북 특성상 모든점이 플랏이 안됨으로 11905개중 랜덤으로 5500개만 선택항여 plot
from folium import IFrame
from folium.plugins import MarkerCluster
import folium
import random as rd
point_num = 5500
m = folium.Map(location=center, tiles="OpenStreetMap", zoom_start=12)
width, height = 310,150
popups, locations = [], []
current_df = candidate_gdf.loc[rd.sample(list(candidate_gdf.index), point_num)]
for idx, row in current_df.iterrows():
locations.append([row['geometry'].y, row['geometry'].x])
h = folium.FeatureGroup(name='candidate')
h.add_children(MarkerCluster(locations=locations))
m.add_children(h)
공공성을 고려하여 현재 배치되지 않은 지역을 기반으로 후보군을 추출하였다.
# 미배치 지역 우선 선정 함수
def join2gid(gdf1, gdf2, target_col, col_name):
"""
input : (candidate_gdf, population_dist_gdfM_not_zero_coverage, "val", "GID_POP")
output : ret_df
describe : gdf1의 좌표가 gdf2의 어느 구역에 들어가는지 파악하여 해당 row의 target_col의 값을 받아온다.
"""
ret_df = copy.deepcopy(gdf1)
start = time.time()
num = len(ret_df)
ret_df[col_name] = 0
for idx, row in ret_df.iterrows():
if idx % 1000 == 0:
print("진행률 : {} / {}".format(idx, num))
print("진행시간 : {}".format(time.time() - start))
temp = gdf2[gpd.GeoSeries.contains(gdf2['geometry'], row.loc['geometry'])][target_col]
if len(temp) > 0:
ret_df[col_name][idx] += float(temp)
return ret_df
# 사용하기 전 변수 세팅
candidate_gdf = candidate_gdf.to_crs(meter_crs)
args = [candidate_gdf, population_dist_gdfM_not_zero_coverage, "val", "GID_POP"]
# 함수 실행
candidate_gdf = join2gid(*args)
candidate_gdf.head()
# 상업지역에 포함되는 100m*100m 지역 뽑아내기
house_commerical_area = house_counts_gdf[house_counts_gdf["commerical_area"] > 0]
house_commerical_area.head()
# 후보군 중에서 상업지역에 들어가는 후보군 파악한다
# 사용하기 전 변수 세팅
args = [candidate_gdf, house_commerical_area, "commerical_area", "commerical_area"]
# 함수 실행
candidate_gdf_commerical = join2gid(*args)
candidate_gdf_commerical.head()
# 후보군 중에서 상업지역에 들어가는 후보군 선별한다.(n=203)
candidate_gdf_commerical_temp = candidate_gdf_commerical[candidate_gdf_commerical["commerical_area"] > 0]
candidate_gdf_commerical_temp.reset_index(drop=True, inplace = True)
del candidate_gdf_commerical_temp["commerical_area"]
candidate_gdf_commerical_temp
# 미배치 지역 중 너무 적은 인구가 거주할 경우(미배치 지역 주거인구 = 200) 제거한다.
gid_pop = 200
candidate_not_cover_df = candidate_gdf[candidate_gdf['GID_POP'] > gid_pop]
candidate_not_cover_df
# 앞서 필터링한 미배치 지역 후보군에 상업지역에 들어가는 후보군 붙여준다
candidate_gdf_with_commerical_area = pd.concat([candidate_not_cover_df, candidate_gdf_commerical_temp], join="outer", axis= 0)
candidate_gdf_with_commerical_area
# 현재까지 제작된 후보군에 200m안 거주 건물 수의 평균을 더해준다.
# 사용하기 전 변수 세팅
candidate_gdf = candidate_gdf.to_crs(meter_crs)
args = [candidate_gdf, house_counts_cover_gdf, 200, gpd.GeoSeries.within, np.sum]
# 함수 실행
candidate_gdf = join_info2station(*args)
candidate_gdf.head() # 결과 확인
# 평균 계산이 되지 않은 경우 nan으로 생성된다. 이를 0으로 변경하자.
candidate_gdf.fillna(0, inplace=True)
candidate_gdf.head()
# 주거 지역 (house_count_mean), 빌딩 연면적, 거주인구 일정 이상의 필터를 적용한다.
pop_count = 1000
b_area_count = 20000
house_count = 10
candidate_pop_gdf = candidate_gdf[
(candidate_gdf['pop_200_sum'] > pop_count) |
(candidate_gdf['b_area_200_sum'] > b_area_count) &
(candidate_gdf['HOUSE_COUNTS_200_mean'] > house_count)
]
candidate_pop_gdf.head()
filter_list = [
candidate_gdf_with_commerical_area,
candidate_pop_gdf,
]
filter_index = []
for df in filter_list:
filter_index += list(df.index)
filter_index = list(set(filter_index))
len(filter_index)
candidate_gdf = candidate_gdf.loc[filter_index]
candidate_gdf.reset_index(drop=True, inplace = True)
candidate_gdf.head()
접근성은 아래와 같이 총 7개의 filter를 적용했다.
candidate_gdf = candidate_gdf.to_crs(default_crs)
candidate_gdf['LAT'] = candidate_gdf['geometry'].apply(lambda x : x.coords[0][1])
candidate_gdf['LON'] = candidate_gdf['geometry'].apply(lambda x : x.coords[0][0])
candidate_gdf.head()
candidate_gdf = countCloseLandmarks(candidate_gdf, busstop_df, 0.3, "busstop")
candidate_gdf.head()
sns.distplot(candidate_gdf['busstop_cnt_300m'])
300m 기준으로 분포한 버스정류장의 개수를 그려보았다.
이 때, 5개가 가장 맡은 값을 보이고 있다. 이점을 기반으로 주변 정류장의 개수가 5개 이상인 후보군을 선택하였다.
bus_count = 5
candidate_bus_count_gdf = candidate_gdf[candidate_gdf['busstop_cnt_300m'] > bus_count]
candidate_gdf = calculateCloseLandmarkFeatures(candidate_gdf, busstop_df, 0.3, "busstop", "GETON_CNT", np.mean)
candidate_gdf.head()
candidate_gdf.fillna(0, inplace=True)
candidate_gdf.head()
sns.distplot(candidate_gdf['busstop_GETON_CNT_mean_300m'])
300m 기준으로 버스 정류장 승하차 인원을 살펴보았을 때, 위와 같은 분포를 보인다.
우리는 승하차 인원이 많은 정류장을 기반으로 후보군을 뽑기 위해 5000명 이상의 승하차 인원을 보이는 후보군을 filter로 적용했다.
bus_mean = 5000
candidate_bus_mean_gdf = candidate_gdf[candidate_gdf['busstop_GETON_CNT_mean_300m'] > bus_mean]
candidate_bus_mean_gdf.head()
candidate_gdf = countCloseLandmarks(candidate_gdf, subway_space_df, 1.0, "subway")
candidate_gdf
sns.distplot(candidate_gdf['subway_cnt_1000m'])
반경 1km 기준 지하철역 개수의 분포는 위와 같다.
우리는 주변에 지하철 역의 유무에 따른 filter를 적용하였다.
subway_cnt = 0
candidate_subway_count_gdf = candidate_gdf[candidate_gdf['subway_cnt_1000m'] > subway_cnt]
candidate_subway_count_gdf.head()
candidate_gdf = calculateCloseLandmarkFeatures(candidate_gdf, subway_space_df, 1.0, "subway", "IN", sum)
candidate_gdf = calculateCloseLandmarkFeatures(candidate_gdf, subway_space_df, 1.0, "subway", "OUT", sum)
candidate_gdf.head()
candidate_gdf['subway_IN+OUT_sum_1000m'] = candidate_gdf['subway_IN_sum_1000m'] + candidate_gdf['subway_OUT_sum_1000m']
candidate_gdf.head()
sns.distplot(candidate_gdf['subway_IN+OUT_sum_1000m'])
주변 1km내 지하철역의 승하차 인원의 합의 분포는 위와 같다.
지하철 역이 있는 경우가 드물기 때문에, 우리는 주변에 지하철역이 존재하기만 한다면 이 승하차인원을 사용하도록 하는 filter를 적용했다.
subway_mean = 0
candidate_subway_mean_gdf = candidate_gdf[candidate_gdf['subway_IN+OUT_sum_1000m'] > subway_mean]
candidate_subway_mean_gdf.head()
candidate_gdf = countCloseLandmarks(candidate_gdf, school_gdf, 1.0, "school")
candidate_gdf.head()
sns.distplot(candidate_gdf['school_cnt_1000m'])
주변 1km내 학교의 개수의 분포는 위와 같다.
우리는 학교가 많은 경우의 후보군을 추출하고 싶었기에 5개 이상의 학교를 가지는 후보군을 추출하는 filter를 적용하였다.
school_count = 5
candidate_school_count_gdf = candidate_gdf[candidate_gdf['school_cnt_1000m'] > school_count]
candidate_school_count_gdf.head()
candidate_gdf = countCloseLandmarks(candidate_gdf, sport_complex_info_df, 1.0, "sports")
candidate_gdf.head()
sns.distplot(candidate_gdf['sports_cnt_1000m'])
주변 1km 체육 시설의 분포는 위와 같다.
우리는 체육시설을 대표하는 후보군을 선정하고 싶었기 때문에 5개 이상의 체육시설을 가지는 후보군을 추출하였다.
sport_count = 5
candidate_sport_count_gdf = candidate_gdf[candidate_gdf['sports_cnt_1000m'] > sport_count]
candidate_sport_count_gdf.head()
candidate_gdf = candidate_gdf.to_crs(default_crs)
candidate_gdf = join_height2station(candidate_gdf, height_gdf)
candidate_gdf.head()
sns.distplot(candidate_gdf['altitude'])
고도에 따른 분포는 위와 같다.
20m 근처에 가장 많은 후보군이 위치하고, 낮은 지역에서 많은 수요가 발생함을 EDA로 부터 확인했다.
이 근거를 바탕으로 22m 보다 작은 후보군을 추출하였다.
min_height = 15
max_height = 30
candidate_altitude_gdf = candidate_gdf[candidate_gdf['altitude'] < max_height]
candidate_altitude_gdf = candidate_altitude_gdf[candidate_altitude_gdf['altitude'] > min_height]
candidate_altitude_gdf
# 위에서 추출한 모든 필터를 합친다.
filter_list = [
candidate_bus_count_gdf,
candidate_bus_mean_gdf,
candidate_subway_count_gdf,
candidate_subway_mean_gdf,
candidate_school_count_gdf,
candidate_sport_count_gdf,
candidate_altitude_gdf
]
filter_index = []
for df in filter_list:
filter_index += list(df.index)
filter_index = list(set(filter_index))
len(filter_index)
candidate_after_filter_gdf = candidate_gdf.loc[filter_index]
candidate_after_filter_gdf
우리는 굉장히 추상적인 후보군 추출을 앞선 EDA의 결과를 종합하여 상당히 가능성이 높은 후보군을 다수 추출하였다. 해당 후보군을 추출하는데 고려한 점은 다음과 같다.
제안된 후보군의 수요를 예측하기 위한 모델을 제작한다.
자전거 공유 서비스를 사용하는데 있어서 수요는 다양한 측면에서 바라볼 수 있다.
이 때, 우리는 이 수요라는 것을 총 두가지 관점에서 바라보았다.
첫번째, 거치대 반입량, 반출량의 절대적인 량을 사용하였다. 이 지표는 하나의 거치대에서 방문하는 사용자의 총량을 대변한다.
두번째, 거치대 반입량, 반출량의 차이를 사용하였다. 이 지표는 하나의 거치대에서 유입, 유출되는 자전거의 변화량을 의미한다.
이 때, 우리는 이러한 두가지 지표의 시간 단위를 일로 설정하였다.
결과적으로 사용하는 수요의 개념은 다음과 같다.
용어 정의
총 자전거 이용량 (반입량+반출량)의 일별 평균 : 절대 수요
총 자전거 변화량 (반입량-반출량)의 일별 평균 : 상대 수요
EDA결과를 기반으로 다양한 변수를 생성하여 사용하였다.
train_df = station_df[['STATION_ID', 'geometry', 'LAT', 'LON']]
train_df.head()
# 거치대 df에 사용할 수 있는 모든 변수를 추가한다.
def createVariableInTrain(train_df):
col_name = ['busstop_cnt_200m', 'busstop_cnt_300m', 'busstop_cnt_500m']
param = [0.2, 0.3, 0.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = countCloseLandmarks(train_df, busstop_df, param[i], "busstop")
col_name = ['busstop_GETON_CNT_sum_200m', 'busstop_GETON_CNT_sum_300m', 'busstop_GETON_CNT_sum_500m']
param = [0.2, 0.3, 0.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = calculateCloseLandmarkFeatures(train_df, busstop_df, param[i], "busstop", "GETON_CNT", sum)
col_name = ['busstop_GETON_CNT_mean_200m', 'busstop_GETON_CNT_mean_300m', 'busstop_GETON_CNT_mean_500m']
param = [0.2, 0.3, 0.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = calculateCloseLandmarkFeatures(train_df, busstop_df, param[i], "busstop", "GETON_CNT", np.mean)
col_name = ['subway_cnt_500m', 'subway_cnt_1000m', 'subway_cnt_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = countCloseLandmarks(train_df, subway_space_df, param[i], "subway")
col_name = ['subway_IN_sum_500m', 'subway_IN_sum_1000m', 'subway_IN_sum_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = calculateCloseLandmarkFeatures(train_df, subway_space_df, param[i], "subway", "IN", sum)
col_name = ['subway_IN_mean_500m', 'subway_IN_mean_1000m', 'subway_IN_mean_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = calculateCloseLandmarkFeatures(train_df, subway_space_df, param[i], "subway", "IN", np.mean)
col_name = ['subway_OUT_sum_500m', 'subway_OUT_sum_1000m', 'subway_OUT_sum_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = calculateCloseLandmarkFeatures(train_df, subway_space_df, param[i], "subway", "OUT", sum)
col_name = ['subway_OUT_mean_500m', 'subway_OUT_mean_1000m', 'subway_OUT_mean_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = calculateCloseLandmarkFeatures(train_df, subway_space_df, param[i], "subway", "OUT", np.mean)
col_name = ['subway_NET_TOTAL_sum_500m', 'subway_NET_TOTAL_sum_1000m', 'subway_NET_TOTAL_sum_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = calculateCloseLandmarkFeatures(train_df, subway_space_df, param[i], "subway", "NET_TOTAL", sum)
col_name = ['subway_NET_TOTAL_mean_500m', 'subway_NET_TOTAL_mean_1000m', 'subway_NET_TOTAL_mean_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = calculateCloseLandmarkFeatures(train_df, subway_space_df, param[i], "subway", "NET_TOTAL", np.mean)
col_name = ['school_cnt_500m', 'school_cnt_1000m', 'school_cnt_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = countCloseLandmarks(train_df, school_gdf, param[i], "school")
col_name = ['sports_cnt_500m', 'sports_cnt_1000m', 'sports_cnt_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(train_df.columns):
train_df = countCloseLandmarks(train_df, sport_complex_info_df, param[i], "sports")
col_name = 'altitude'
print("{} 작업중".format(col_name))
if col_name not in list(train_df.columns):
train_df = join_height2station(train_df, height_gdf)
train_df = train_df.to_crs(meter_crs)
col_name = 'GID_POP'
print("{} 작업중".format(col_name))
if col_name not in list(train_df.columns):
# 사용하기 전 변수 세팅
args = [train_df, population_dist_gdfM_not_zero_coverage, "val", "GID_POP"]
# 함수 실행
train_df = join2gid(*args)
col_name = 'pop_200_sum'
print("{} 작업중".format(col_name))
if col_name not in list(train_df.columns):
# 사용하기 전 변수 세팅
args = [train_df, house_counts_cover_gdf, 200, gpd.GeoSeries.within, np.sum]
# 함수 실행
train_df = join_info2station(*args)
col_name = 'areatype'
print("{} 작업중".format(col_name))
if col_name not in list(train_df.columns):
# 사용하기 전 변수 세팅
args = [train_df, house_counts_gdf, "areatype", "areatype"]
# 함수 실행
train_df = join2gid(*args)
train_df.fillna(0, inplace=True)
col_name = 'CLUSTER'
print("{} 작업중".format(col_name))
if col_name not in list(train_df.columns):
cluster = pd.DataFrame([list(partition.keys()), list(partition.values())]).T
cluster.columns = ['STATION_ID', 'CLUSTER']
train_df = pd.merge(train_df, cluster, how='left', on='STATION_ID')
train_df.dropna(inplace=True)
train_df = train_df.to_crs(default_crs)
return train_df
train_df = createVariableInTrain(train_df)
train_df
# 수요로 사용할 수 있는 변수 가져오기
station_demand_df.columns = ['STATION_ID',
'LEAS_COUNTS',
'RTN_COUNTS',
'TOTAL',
'RTN-LEAS',
'LEAS_MEAN',
'RTN_MEAN',
'TOTAL_MEAN',
'RTN-LEAS_MEAN',
]
station_demand_df
station_demand_df.fillna(0, inplace=True)
# train df에 이 이용량 수요 데이터를 추가해준다.
train_df = pd.merge(train_df, station_demand_df, how='left', on='STATION_ID')
train_df.dropna(inplace=True)
train_df.head()
# train 데이터 info
# 훈련에 앞서 실수형 데이터 아닌 변수가 있는지 확인
train_df.info()
# 훈련에 필요한 컬럼을 찾는다.
remove_list = ['STATION_ID','geometry','LAT','LON',
'LEAS_COUNTS', 'RTN_COUNTS', 'TOTAL', 'RTN-LEAS', 'LEAS_MEAN', 'RTN_MEAN', 'TOTAL_MEAN', 'RTN-LEAS_MEAN']
data_columns = list(train_df.columns)
for rmv in remove_list:
data_columns.remove(rmv)
data_columns
extratree, lightgbm, xgboost로 를 테스트 해보았고, 그중 가장 좋은 성능을 보인 extratree를 사용하였다. 훈련 시간 관계로 extratree만 훈련하도록 하였다.
# 0번 스테이션의 평균값을 더해준다.
train_df["RTN-LEAS_MEAN"] += bias
model_out = dict()
import xgboost
import lightgbm
MLA = [
ensemble.ExtraTreesRegressor()
#lightgbm.LGBMRegressor(),
#xgboost.XGBRegressor()
]
MLA
Targets = ['TOTAL_MEAN', 'RTN-LEAS_MEAN']
for Target in Targets:
data = copy.deepcopy(train_df)
alg = copy.deepcopy(MLA[0])
alg.fit(data[data_columns], data[Target])
model_out[Target] = copy.deepcopy(alg)
pred = alg.predict(data[data_columns])
print('-------------------------')
print(Target , 'loss', ((pred - data[Target].values)**2).sum())
print('-------------------------')
현재 설치된 거치대를 기준으로 절대수요, 상대수요를 예측한 결과 높은 성능의 모델이 도출되었다. 데이터 개수가 많이 부족하지만, 상대적으로 가벼운 모델인 decision tree 기반 모델을 사용하여 최대한 일반화 성능을 높히는 방향으로 모델링을 진행하였다.
이제 만들어진 후보군을 기반으로 훈련한 모델에 입력함으로써 절대수요와 상대수요를 예측한다.
앞서 진행한 EDA에서 네트워크 분석을 통해 거치대를 총 6개의 군집으로 나누었다. 우리가 제작한 후보군을 해당 군집에 적용한다. 이 때, KNN 알고리즘을 사용하여 유클리드 거리 기반으로 가까운 군집 번호를 부여한다.
4.3.2에서 이 클러스터를 배정해 주기 위해, 기존 거치대의 정보를 KNN 알고리즘에 설정한다.

train_knn = station_df.dropna()[['geometry', 'CLUSTER']]
train_knn = train_knn.to_crs(meter_crs)
train_knn
# KNN을 훈련하기 위한 훈련 데이터를 제작한다.
training_x1 = np.array(train_knn['geometry'].apply(lambda x : x.coords[0][0])).reshape(-1, 1)
training_x2 = np.array(train_knn['geometry'].apply(lambda x : x.coords[0][1])).reshape(-1, 1)
training_x = np.concatenate((training_x1, training_x2), axis=1)
training_x
# 훈련 데이터의 군집 번호를 뽑아낸다.
training_y = np.array(train_knn['CLUSTER'].apply(lambda x : int(x)))
training_y
# 위의 두 데이터를 기반으로 KNN 모델을 제작해 둔다.
# 이제 특정 후보군이 들어왔을 때, 가장 가까운 점 3개를 기반으로 네트워크 군집을 배정한다.
from sklearn.neighbors import KNeighborsClassifier
KNN = KNeighborsClassifier(n_neighbors = 3)
KNN.fit(training_x, training_y)
수요 예측 모델을 훈련했을 때 사용한 모든 변수를 후보군에 붙인다.
candidate_after_filter_gdf.reset_index(drop=True, inplace=True)
test_df = copy.deepcopy(candidate_after_filter_gdf)
test_df
def createVariableInTest(test_df):
col_name = ['busstop_cnt_200m', 'busstop_cnt_300m', 'busstop_cnt_500m']
param = [0.2, 0.3, 0.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = countCloseLandmarks(test_df, busstop_df, param[i], "busstop")
col_name = ['busstop_GETON_CNT_sum_200m', 'busstop_GETON_CNT_sum_300m', 'busstop_GETON_CNT_sum_500m']
param = [0.2, 0.3, 0.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = calculateCloseLandmarkFeatures(test_df, busstop_df, param[i], "busstop", "GETON_CNT", sum)
col_name = ['busstop_GETON_CNT_mean_200m', 'busstop_GETON_CNT_mean_300m', 'busstop_GETON_CNT_mean_500m']
param = [0.2, 0.3, 0.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = calculateCloseLandmarkFeatures(test_df, busstop_df, param[i], "busstop", "GETON_CNT", np.mean)
col_name = ['subway_cnt_500m', 'subway_cnt_1000m', 'subway_cnt_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = countCloseLandmarks(test_df, subway_space_df, param[i], "subway")
col_name = ['subway_IN_sum_500m', 'subway_IN_sum_1000m', 'subway_IN_sum_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = calculateCloseLandmarkFeatures(test_df, subway_space_df, param[i], "subway", "IN", sum)
col_name = ['subway_IN_mean_500m', 'subway_IN_mean_1000m', 'subway_IN_mean_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = calculateCloseLandmarkFeatures(test_df, subway_space_df, param[i], "subway", "IN", np.mean)
col_name = ['subway_OUT_sum_500m', 'subway_OUT_sum_1000m', 'subway_OUT_sum_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = calculateCloseLandmarkFeatures(test_df, subway_space_df, param[i], "subway", "OUT", sum)
col_name = ['subway_OUT_mean_500m', 'subway_OUT_mean_1000m', 'subway_OUT_mean_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = calculateCloseLandmarkFeatures(test_df, subway_space_df, param[i], "subway", "OUT", np.mean)
col_name = ['subway_NET_TOTAL_sum_500m', 'subway_NET_TOTAL_sum_1000m', 'subway_NET_TOTAL_sum_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = calculateCloseLandmarkFeatures(test_df, subway_space_df, param[i], "subway", "NET_TOTAL", sum)
col_name = ['subway_NET_TOTAL_mean_500m', 'subway_NET_TOTAL_mean_1000m', 'subway_NET_TOTAL_mean_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = calculateCloseLandmarkFeatures(test_df, subway_space_df, param[i], "subway", "NET_TOTAL", np.mean)
col_name = ['school_cnt_500m', 'school_cnt_1000m', 'school_cnt_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = countCloseLandmarks(test_df, school_gdf, param[i], "school")
col_name = ['sports_cnt_500m', 'sports_cnt_1000m', 'sports_cnt_1500m']
param = [0.5, 1.0, 1.5]
for i in range(len(col_name)):
print("{} 작업중".format(col_name[i]))
if col_name[i] not in list(test_df.columns):
test_df = countCloseLandmarks(test_df, sport_complex_info_df, param[i], "sports")
col_name = 'altitude'
print("{} 작업중".format(col_name))
if col_name not in list(test_df.columns):
test_df = join_height2station(test_df, height_gdf)
test_df = test_df.to_crs(meter_crs)
col_name = 'GID_POP'
print("{} 작업중".format(col_name))
if col_name not in list(test_df.columns):
# 사용하기 전 변수 세팅
args = [test_df, population_dist_gdfM_not_zero_coverage, "val", "GID_POP"]
# 함수 실행
test_df = join2gid(*args)
col_name = 'pop_200_sum'
print("{} 작업중".format(col_name))
if col_name not in list(test_df.columns):
# 사용하기 전 변수 세팅
args = [test_df, house_counts_cover_gdf, 200, gpd.GeoSeries.within, np.sum]
# 함수 실행
test_df = join_info2station(*args)
col_name = 'areatype'
print("{} 작업중".format(col_name))
if col_name not in list(test_df.columns):
# 사용하기 전 변수 세팅
args = [test_df, house_counts_gdf, "areatype", "areatype"]
# 함수 실행
test_df = join2gid(*args)
col_name = 'CLUSTER'
print("{} 작업중".format(col_name))
if col_name not in list(test_df.columns):
test_x1 = np.array(test_df['geometry'].apply(lambda x : x.coords[0][0])).reshape(-1, 1)
test_x2 = np.array(test_df['geometry'].apply(lambda x : x.coords[0][1])).reshape(-1, 1)
test_x = np.concatenate((test_x1, test_x2), axis=1)
test_df['CLUSTER'] = KNN.predict(test_x)
test_df.fillna(0, inplace=True)
test_df = test_df.to_crs(default_crs)
return test_df
test_df = createVariableInTest(test_df)
test_df
test_data = test_df[data_columns]
test_data
후보군을 기반으로 예상되는 절대 수요와 상대 수요를 구한다.
for key in model_out.keys():
alg = model_out[key]
pred = alg.predict(test_data)
test_df[key] = pred
test_df[['geometry', 'TOTAL_MEAN', 'RTN-LEAS_MEAN']]
예측된 절대 수요와 상대 수요를 기반으로 수요 증대와 Rebalancing을 고려한 최적 후보군을 선정한다.

해당 알고리즘을 구현하는 데 있어 고민한 점은 크게 4가지이다.
가정 : 자전거의 이동 범위는 지역을 크게 벗어나지 않는다.
지역별 요구 거치대 비율을 구현하는 데 있어 고민한 점은 크게 2가지이다.
위의 두가지 요인을 모두 고려한 지역별 요구 거치대 비율 이라는 변수를 생성하였다. 위 변수를 제작한 수식은 다음과 같다.
$normalization({1\over{stationCount}} \times populationRatio)$
station_count가 크면, 거치대는 적게 설치해야 한다. 따라서 역수 처리 하였다. population_ratio가 크면, 거치대는 더 설치해야 한다. 따라서 그대로 두었다.
두 지표를 곱한 후에 이를 정규화하여 해당 지역의 요구 거치대 비율로 사용하였다.
gdf_boundary_town.columns
# 1. station_count가 작으면 값을 크게 -> 역수 예외처리 0일 때는 0.4로 준다.
# 너무 많은 값을 주면 영향력이 커진다.
# 2. total pop이 크면 우선 배정한다. -> 확률로 만들어 둔다.
# 3. 두개를 곱한다.
# 4. 정규화 한다.
total_population = sum(gdf_boundary_town["TOTAL_POP"])
gdf_boundary_town["PROB_STATION_COUNT"] = gdf_boundary_town["STATION_COUNT"].apply(lambda x: np.float(1/x) if x != 0 else 0.4)
gdf_boundary_town["PROB_TOTAL_POP"] = gdf_boundary_town["TOTAL_POP"].apply(lambda x: x/total_population)
gdf_boundary_town["LOCATE_PROB"] = gdf_boundary_town["PROB_STATION_COUNT"] * gdf_boundary_town["PROB_TOTAL_POP"]# * gdf_boundary_town["PROB_STATION_TOTAL_MEAN"]
total_prob = sum(gdf_boundary_town["LOCATE_PROB"])
gdf_boundary_town["LOCATE_PROB"] = gdf_boundary_town["LOCATE_PROB"].apply(lambda x: x/total_prob)
gdf_boundary_town["LOCATE_NUM"] = gdf_boundary_town["LOCATE_PROB"].apply(lambda x: round(x*140))
print("배치 요구 개수 : ", sum(gdf_boundary_town["LOCATE_NUM"]))
gdf_boundary_town[["DONG_NM", "LOCATE_PROB", "LOCATE_NUM","STATION_COUNT", ]]
최종적으로 예측된 절대 수요와 상대 수요를 기반으로 수요 증대와 Rebalancing을 해결하는 장소를 찾는다.
5.1에서 제작한 지역별 요구 거치대 비율을 기반으로 진행하기에 앞서, 이를 쉽게 반영하기 위한 데이터 프레임을 제작한다.
| variable | means |
|---|---|
| total | 최종 행정동에 배치될 거치대의 총 개수 |
| exist | 행정동 내에 존재하는 거치대의 개수 |
| need | 요구 거치대 비율을 기반으로 필요한 거치대의 개수 |
| min_x | 해당 지역의 x 하한 |
| max_x | 해당 지역의 x 상한 |
| min_y | 해당 지역의 y 하한 |
| max_y | 해당 지역의 y 상한 |
x, y에 대한 정보는 추후 Metric을 구하는데 사용될 것이다.
import time
def join_Town2Candidate(gdf1, gdf2, target_col):
"""
input : (candidate_gdf, gdf_boundary_town, "DONG_NM")
output : ret_df
describe : gdf1의 좌표가 gdf2의 행정동을 받아온다.
"""
ret_df = copy.deepcopy(gdf1)
start = time.time()
num = len(ret_df)
ret_df[target_col] = 0
for idx, row in ret_df.iterrows():
if idx % 1000 == 0:
print("진행률 : {} / {}".format(idx, num))
print("진행시간 : {}".format(time.time() - start))
temp = gdf2[gpd.GeoSeries.contains(gdf2['geometry'], row.loc['geometry'])][target_col]
if len(temp) > 0:
ret_df[target_col][idx] = list(temp)[0]
return ret_df
station_dong_df = train_df[['STATION_ID', 'geometry', 'RTN-LEAS_MEAN']]
station_dong_df
# 사용하기 전 변수 세팅
station_dong_df = station_dong_df.to_crs(meter_crs)
args = [station_dong_df, gdf_boundary_town, "DONG_NM"]
# 함수 실행
station_dong_df = join_Town2Candidate(*args)
station_dong_df.head()
# 모든 동이름을 index로 사용한 df를 만든다.
dong_name = list(gdf_boundary_town['DONG_NM'].unique())
dong_df = gpd.GeoDataFrame(columns=['total', 'exist', 'need', 'min_x', 'max_x', 'min_y', 'max_y', 'geometry'], index=dong_name).set_crs(meter_crs)
dong_df.head()
# 각각의 geometry에서 상한 하한을 구해 넣어준다.
for idx, row in gdf_boundary_town.iterrows():
x, y = row['geometry'][0].exterior.xy
min_x = min(x)
max_x = max(x)
min_y = min(y)
max_y = max(y)
dong_df.loc[row['DONG_NM']] = [0, 0, 0, str(min_x), str(max_x), str(min_y), str(max_y), row['geometry']]
dong_df['min_x'] = dong_df['min_x'].apply(lambda x : float(x))
dong_df['max_x'] = dong_df['max_x'].apply(lambda x : float(x))
dong_df['min_y'] = dong_df['min_y'].apply(lambda x : float(x))
dong_df['max_y'] = dong_df['max_y'].apply(lambda x : float(x))
dong_df.head()
# 해당 행정동 내에 현재 존재하는 거치대의 개수를 넣어준다.
dong_df['exist'] = station_dong_df[['DONG_NM', 'STATION_ID']].groupby(['DONG_NM']).count()
dong_df.fillna(0, inplace=True)
dong_df.head()
# 지역별 요구 거치대 비율을 기반으로 필요한 개수를 추가한다.
dong_df["need"] = list(gdf_boundary_town["LOCATE_NUM"])
dong_df["total"] = dong_df["exist"] + dong_df["need"]
dong_df
# 사용하기 전 변수 세팅
test_df = test_df.to_crs(meter_crs)
args = [test_df, gdf_boundary_town, "DONG_NM"]
# 함수 실행
test_df = join_Town2Candidate(*args)
test_df.head()
절대 수요를 기반으로 필요한 후보군의 개수 x 군집 변수(3)에 해당하는 개수의 클러스터링을 진행한다. 지역 별로 요구되는 거치대 개수x3에 해당하는 군집으로 후보군을 분리함으로써, 비슷한 절대 수요를 가지는 후보군을 하나로 묶어, 해당 지역의 특징을 반영한 대표군을 만들고자 하였다.
상대 수요는 간단하게 말하면, 해당 거치대에 반입량의 순수합이다. 즉 양수인 경우 반입이 많고, 음수인 경우 반출이 많다. 거치대의 이러한 특징은 시스템 운영 비용에 직접적으로 영향을 미친다. 따라서 이러한 상대 수요의 총 합이 0이 되는, 즉 Rebalancing에 도움을 주는 지역을 선정하는 것이 매우 중요하다.
위에서 제작한 대표군은 최종적으로 배치될 지역 후보군이다. 하지만 공유 자전거 문제의 특징상 하나의 새 거치대 배치는 지역내에 존재하는 다른 모든 자전거에 영향을 미친다. 이러한 의존성이 있는 문제에 대한 해답으로 우리는 상대 수요를 기반으로 한 이변량 정규분포의 부피로 제안한 거치대의 Rebalancing 영향도를 측정하였다.
$mertic = \int_{y_{min}}^{y_{max}} \int_{x_{min}}^{x_{max}} \left |\Phi ( x_{1}, x_{2})\right |$
$\Phi(x_{1}, x_{2}) = \sum_{i}^{n}\alpha_i *(\frac{exp^{-\frac{1}{2}(x-\mu_{i})^{t}\Sigma_{i}(x-\mu_{i})}}{\sqrt{2\pi \left | \Sigma_{i} \right |}})$
여기서 $\alpha$는 현재 제안된 후보군을 포함한 상대 수요의 총합이다. 이렇게 정의한다면, 이 Metric 함수(= 상대 수요의 총합 x 이변량 정규 분포)의 부피는 지역 전역에 있는 상대 수요의 총량을 대변한다고 할 수 있다. 이러한 지표를 사용하였을 때, 상대 수요의 총량이 0이 되도록 하는 대표 후보군을 선택하는 방법으로 최적 거치대를 선정하였다. 해당 알고리즘에서는 시간 복잡도 문제로 Greedy알고리즘을 사용하였다.
# Set 지역, 지역 비율, 이미 배치된 수, 총 배치대 수
# Set 후보군, 지역 후보군, 지역 후보군 절대 수요
# Set 배치된 거치대
# Set 최종 후보군
# for 지역:
# 추가 배치대 수 = 지역 비율 * 총 배치대 수 - 이미 배치된 수
# 클러스터 수 = 3 * 추가 배치대 수
# 클러스터 = Kmeans(cluster = 클러스터수, weight = 지역 후보군 절대 수요)
# 지역 후보군 = 클러스터.중심
# for 추가 배치대 수:
# 최소 리벨런싱 = INF
# for 지역 후보군:
# 배치된 거치대.push_back(지역 후보군[i])
# 최소 리벨런싱 = MIN(최소 리벨런싱, Metric(배치된 거치대))
# 배치된 거치대.pop_back(지역 후보군[i])
# 배치된 거치대.push_back(지역 후보군[min_index])
# 지역 후보군.remove(지역 후보군[min_index])
# END
# END
# 최종 후보군.push_back(배치된 거치대)
# END
# 실제 작동 코드
# 클러스터링 하기전 모든 후보군
for idx, row in dong_df.iterrows():
exist_station = station_dong_df[station_dong_df['DONG_NM'] == idx]
for i, r in exist_station.iterrows():
x, y = r['geometry'].xy
plt.plot(x, y, 'ro')
candidate_station = test_df[test_df['DONG_NM'] == idx]
for i, r in candidate_station.iterrows():
x, y = r['geometry'].xy
plt.plot(x, y, 'bo')
x, y = row['geometry'][0].exterior.xy
plt.plot(x, y)
클러스터링 하기전의 후보군과 이미 배치된 거치대를 시각화해보았다.
| 색 | 의미 |
|---|---|
| 파랑 | 배치될 수 있는 후보군 |
| 빨강 | 배치된 거치대 |
# 절대 수요 기반 클러스터링 진행 후 압축된 후보군 시각화
def createClusterDf(candidate_station, need):
num = need * 3
if len(candidate_station) < num:
candidate_station['X'] = candidate_station['geometry'].apply(lambda x : x.coords[0][0])
candidate_station['Y'] = candidate_station['geometry'].apply(lambda x : x.coords[0][1])
return candidate_station[['X', 'Y', 'geometry']]
kmeans = sklearn.cluster.KMeans(n_clusters=num)
kmeans.fit(
X=np.stack([candidate_station['geometry'].apply(lambda x : x.coords[0][0]),
candidate_station['geometry'].apply(lambda x : x.coords[0][1])], axis=1),
sample_weight=candidate_station['TOTAL_MEAN'],
)
cluster_df = kmeans.cluster_centers_
cluster_df = gpd.GeoDataFrame(kmeans.cluster_centers_, columns=['X', 'Y'])
cluster_df['XY'] = cluster_df['X'].apply(lambda x : str(x)) + ' ' + cluster_df['Y'].apply(lambda x : str(x))
cluster_df['geometry'] = cluster_df['XY'].apply(lambda x : Point(float(x.split(' ')[0]), float(x.split(' ')[1])))
cluster_df = cluster_df.set_crs(meter_crs)
cluster_df = cluster_df[['X', 'Y', 'geometry']]
return cluster_df
for idx, row in dong_df.iterrows():
exist_station = station_dong_df[station_dong_df['DONG_NM'] == idx]
for i, r in exist_station.iterrows():
x, y = r['geometry'].xy
plt.plot(x, y, 'ro')
if row['need'] != 0:
candidate_station = test_df[test_df['DONG_NM'] == idx]
candidate_station = createClusterDf(candidate_station, row['need'])
for i, r in candidate_station.iterrows():
x, y = r['geometry'].xy
plt.plot(x, y, 'bo')
x, y = row['geometry'][0].exterior.xy
plt.plot(x, y, 'g')
절대 수요 기반 클러스터링을 진행한 후의 압축된 후보군과 배치된 거치대를 비교해보았다.
| 색 | 의미 |
|---|---|
| 파랑 | 배치될 수 있는 압축된 후보군 |
| 빨강 | 배치된 거치대 |
결과적으로 절대 수요 기반 클러스터링의 결과가 기존의 후보군을 잘 대표해주며, 이를 압축한 결과를 보여주었다. 이를 기반으로 해당 알고리즘을 수행하게 된다.
# 고양시 상하한 설정 -> 추후 부피를 구하기 위함
min_x = np.inf; max_x = 0;
min_y = np.inf; max_y = 0;
for idx, row in dong_df.iterrows():
x, y = row['geometry'][0].exterior.xy
min_x = min(min_x , min(x))
max_x = max(max_x , max(x))
min_y = min(min_y , min(y))
max_y = max(max_y , max(y))
print("min_x : ", min_x)
print("max_x : ", max_x)
print("min_y : ", min_y)
print("max_y : ", max_y)
# Rebalancing을 고려한 최적 후보군 선정 알고리즘
import scipy.integrate
from scipy.integrate import dblquad
from numpy import exp,arange
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show
def createClusterDf(candidate_station, need):
num = need * 3
if len(candidate_station) < num:
candidate_station['X'] = candidate_station['geometry'].apply(lambda x : x.coords[0][0])
candidate_station['Y'] = candidate_station['geometry'].apply(lambda x : x.coords[0][1])
candidate_station.reset_index(drop=True, inplace=True)
return candidate_station[['X', 'Y', 'geometry']]
kmeans = sklearn.cluster.KMeans(n_clusters=num)
kmeans.fit(
X=np.stack([candidate_station['geometry'].apply(lambda x : x.coords[0][0]),
candidate_station['geometry'].apply(lambda x : x.coords[0][1])], axis=1),
sample_weight=candidate_station['TOTAL_MEAN'],
)
cluster_df = kmeans.cluster_centers_
cluster_df = gpd.GeoDataFrame(kmeans.cluster_centers_, columns=['X', 'Y'])
cluster_df['XY'] = cluster_df['X'].apply(lambda x : str(x)) + ' ' + cluster_df['Y'].apply(lambda x : str(x))
cluster_df['geometry'] = cluster_df['XY'].apply(lambda x : Point(float(x.split(' ')[0]), float(x.split(' ')[1])))
cluster_df = cluster_df.set_crs(meter_crs)
cluster_df = cluster_df[['X', 'Y', 'geometry']]
return cluster_df
def predData(candidate_station, data_columns, model_out):
candidate_station = candidate_station.to_crs(default_crs)
candidate_station['LAT'] = candidate_station['geometry'].apply(lambda x : x.coords[0][1])
candidate_station['LON'] = candidate_station['geometry'].apply(lambda x : x.coords[0][0])
candidate_station = candidate_station[['LON', 'LAT', 'geometry']]
candidate_station = createVariableInTest(candidate_station)
test_data = candidate_station[data_columns]
alg = model_out["RTN-LEAS_MEAN"]
pred = alg.predict(test_data)
candidate_station["RTN-LEAS_MEAN"] = pred
candidate_station = candidate_station.to_crs(meter_crs)
candidate_station['X'] = candidate_station['geometry'].apply(lambda x : x.coords[0][0])
candidate_station['Y'] = candidate_station['geometry'].apply(lambda x : x.coords[0][1])
return candidate_station[['X', 'Y', 'RTN-LEAS_MEAN', 'geometry']]
def f(x, y):
global X, Y
sigma = 1.0
scale = 1
z = 0
for i in range(X.shape[0]):
rv = sp.stats.multivariate_normal([X[i, 0], X[i, 1]], [[sigma * sigma, 0.0], [0.0, sigma * sigma]])
z += scale*(Y[i])*rv.pdf([x, y])
return abs(z)
final_data_x = np.array([[0, 0]])
final_data_y = np.array([0])
num = 0
for idx, row in dong_df.iterrows():
start = time.time()
print("---------------",idx,"---------------")
print("{} / {}".format(num, len(dong_df)))
num += 1
exist_station = station_dong_df[station_dong_df['DONG_NM'] == idx]
candidate_station = test_df[test_df['DONG_NM'] == idx]
if row['need'] == 0:
if len(exist_station['geometry']) > 0:
x1 = np.array(exist_station['geometry'].apply(lambda x : x.coords[0][0])).reshape(-1, 1)
x2 = np.array(exist_station['geometry'].apply(lambda x : x.coords[0][1])).reshape(-1, 1)
X = np.hstack((x1, x2))
Y = np.array(exist_station['RTN-LEAS_MEAN'])
final_data_x = np.vstack((final_data_x, X))
final_data_y = np.hstack((final_data_y, Y))
continue
candidate_station = createClusterDf(candidate_station, row['need'])
if len(candidate_station) == 0:
continue
candidate_station = predData(candidate_station, data_columns, model_out)
candidate_station.sort_values(by=['RTN-LEAS_MEAN'], axis=0, ascending=False, inplace=True)
candidate_station.reset_index(drop=True, inplace=True)
x1 = np.array(exist_station['geometry'].apply(lambda x : x.coords[0][0])).reshape(-1, 1)
x2 = np.array(exist_station['geometry'].apply(lambda x : x.coords[0][1])).reshape(-1, 1)
X = np.hstack((x1, x2))
Y = np.array(exist_station['RTN-LEAS_MEAN'])
ok = []
if len(candidate_station) < row['need']:
row['need'] = len(candidate_station)
for t in range(row['need']):
min_value = np.inf
min_idx = t
param = []
for i, r in candidate_station.iterrows():
print("reblancing 진행중 {} / {} || {} / {}".format(i, len(candidate_station), t, row['need']))
if i in ok:
continue
X = np.vstack((X, np.array([[r['X'], r['Y']]])))
Y = np.hstack((Y, np.array(r['RTN-LEAS_MEAN'])))
#area = dblquad(lambda x, y : f(x, y), min_y, max_y, min_x, max_x)
#area = dblquad(lambda x, y : f(x, y), min_x, max_x, min_y, max_y)
area = dblquad(lambda x, y : f(x, y), row['min_x'], row['max_x'], row['min_y'], row['max_y'])
X = X[:-1, :]
Y = Y[:-1]
if area[0] < min_value:
min_value = area[0]
min_idx = i
param = [r['X'], r['Y'], r['RTN-LEAS_MEAN']]
X = np.vstack((X, np.array([[param[0], param[1]]])))
Y = np.hstack((Y, np.array(param[2])))
ok.append(min_idx)
final_data_x = np.vstack((final_data_x, X))
final_data_y = np.hstack((final_data_y, Y))
print("시간 : ", time.time() - start)
final_data_x = final_data_x[1:, :]
final_data_y = final_data_y[1:]
# Rebalancing을 해결하는 최적 후보군 시각화
for idx, row in dong_df.iterrows():
x, y = row['geometry'][0].exterior.xy
plt.plot(x, y)
# axis("off")
plt.plot(final_data_x[:, 0], final_data_x[:, 1], 'bo')
전체적으로 거치대가 잘 뿌려진것을 확인할수 있다.
# 정규 분포계산 함수
def calc_f(x, y):
global X, Y
sigma = 300.0
scale = 10000
if len(x.shape) == 1:
z = np.zeros((x.shape[0]))
elif len(x.shape) == 2:
z = np.zeros((x.shape[0], x.shape[1]))
pos = np.dstack((x, y))
num = X.shape[0]
for i in range(X.shape[0]):
if i % 50 == 0:
print("진행중 {} / {}".format(i, num))
rv = sp.stats.multivariate_normal([X[i, 0], X[i, 1]], sigma*sigma*np.eye(2))
z += scale*(Y[i])*rv.pdf(pos)
return z
# 고양시 전체를 기준으로 z값 계산
X = final_data_x
Y = final_data_y
x, y = np.mgrid[min_x:max_x:10, min_y:max_y:10]
z = calc_f(x, y)
#고양시 전체 Rebalancing 정도
fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection='3d')
ax.contour3D(x, y, z, 100, cmap=plt.cm.rainbow)
plt.title("Relative Demand 3D Plot", fontsize=20)
plt.show()
극단적으로 치운친 봉오리가 많이 없는 것으로 보아 리벨런싱이 잘 되었다고 생각할 수 있다.
하지만 중간에 극단적인 부분이 있다. 이 부분을 확대해서 관찰할 필요성이 있다.
# 후보군 포함 상대 수요 히트맵
from numpy import exp,arange
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show
# fig, ax = plt.subplots()
fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(1,1,1)
im = imshow(z.T,cmap=plt.cm.rainbow)
cset = ax.contour(z.T,arange(-1,1.5,0.2),linewidths=2,cmap=cm.Set2)
ax.clabel(cset,inline=True,fmt='%1.1f',fontsize=10)
ax.axis('off')
colorbar(im, shrink=0.7)
plt.title("Relative Demand Heatmap In Goyang City", fontsize=20)
중간에 극단적인 부분이 호수공원, 킨텍스, 대화역, 주엽역, 라페스타, 웨스턴돔 부분인 것을 확인했다.
해당 부분만 떼어내서 관찰해보자.
# 호수공원, 킨텍스, 대화역, 주엽역, 라페스타, 웨스턴돔 부분 plot
row_index = ['주엽1동', '주엽2동', '대화동','정발산동','중산동','탄현동', '일산1동', '일산2동', '일산3동']
min_x_ = dong_df.loc[row_index]['min_x'].min()
max_x_ = dong_df.loc[row_index]['max_x'].max()
min_y_ = dong_df.loc[row_index]['min_y'].min()
max_y_ = dong_df.loc[row_index]['max_y'].max()
x, y = np.mgrid[min_x_:max_x_:10, min_y_:max_y_:10]
z = calc_f(x, y)
# 호수공원, 킨텍스, 대화역, 주엽역, 라페스타, 웨스턴돔 부분 plot
fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection='3d')
ax.contour3D(x, y, z, 1000, cmap=plt.cm.rainbow)
plt.title("Relative Demand 3D plot Of Minus Region With Candidates", fontsize=20)
plt.show()
기존의 설치된 거치대의 상대수요와 Rebalancing된 후의 상태를 비교해보자.
x1 = np.array(station_dong_df['geometry'].apply(lambda x : x.coords[0][0])).reshape(-1, 1)
x2 = np.array(station_dong_df['geometry'].apply(lambda x : x.coords[0][1])).reshape(-1, 1)
X = np.hstack((x1, x2))
Y = np.array(station_dong_df['RTN-LEAS_MEAN'] + bias)
row_index = ['주엽1동', '주엽2동', '대화동','정발산동','중산동','탄현동', '일산1동', '일산2동', '일산3동']
min_x_ = dong_df.loc[row_index]['min_x'].min()
max_x_ = dong_df.loc[row_index]['max_x'].max()
min_y_ = dong_df.loc[row_index]['min_y'].min()
max_y_ = dong_df.loc[row_index]['max_y'].max()
x, y = np.mgrid[min_x_:max_x_:10, min_y_:max_y_:10]
z = calc_f(x, y)
# 기존 거치대 기반 상대수요 시각화
fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection='3d')
ax.contour3D(x, y, z, 1000, cmap=plt.cm.rainbow)
plt.title("Relative Demand 3D plot Of Minus Region With Installed Station", fontsize=20)
plt.show()
후보군을 포함하지 않았을 때, 상대 수요의 분포를 보면 굉장히 들쑥날쑥함을 알 수 있다.
제안한 후보군을 추가한 이후 그래프가 완만해졌고, 이는 곧 어느정도 불균형을 해소했음을 의미한다.
하지만 여전히 거치대에 따른 상대 수요의 불균형이 존재하고, 이를 해결하기 위해 자전거 운송 차량을 배치할 필요성이 엿보인다.
마지막으로 제안된 후보군의 위치를 보정해주는 작업을 진행한다.
가장 가까운 도로로 후보군들을 이동시킨다. 그리고 마지막으로 답안 양식으로 후보군을 추가한다.
# 후보군의 geometry 정보만을 가져온다.
final_gdf = gpd.GeoDataFrame(final_data_x, columns=['X', 'Y'])
final_gdf['XY'] = final_gdf['X'].apply(lambda x : str(x)) + ' ' + final_gdf['Y'].apply(lambda x : str(x))
final_gdf['geometry'] = final_gdf['XY'].apply(lambda x : Point(float(x.split()[0]), float(x.split()[1])))
final_gdf = final_gdf.set_crs(meter_crs)
final_gdf = final_gdf.to_crs(default_crs)
final_gdf = final_gdf[['geometry']]
final_gdf
# 도로 데이터에서 geometry 정보를 제외한 describe
address_road_gdf.drop('geometry', axis=1).describe(include='all')
# 모든 도로의 RDS_MAN_NO를 형변환한다.
for column in ['RDS_MAN_NO']:
address_road_gdf[column] = address_road_gdf[column].apply(str)
address_road_gdf.drop('geometry', axis=1).describe(include='all')
# 도로 정보를 점으로 변경한다.
import shapely
x0, x1 = list(), list()
for lines in address_road_gdf['geometry']:
for line in lines:
for _x0, _y0, _x1, _y1 in zip(line.xy[0][:-1], line.xy[1][:-1], line.xy[0][1:], line.xy[1][1:]):
x0.append([_x0, _y0])
x1.append([_x1, _y1])
lines = shapely.ops.cascaded_union([shapely.geometry.LineString([x0, x1]) for x0, x1 in zip(x0, x1)])
# 가장 가까운 점을 가져온다.
def find_nearest_point(x0, x1, x2):
x0 = tf.reshape(tf.constant(x0, dtype=tf.float32), [-1, 2, 1])
x1 = tf.reshape(tf.constant(x1, dtype=tf.float32), [-1, 2, 1])
x2 = tf.reshape(tf.constant(x2, dtype=tf.float32), [-1, 2, 1])
A = tf.reshape(
tf.stack([
(x0[:, 1, 0] - x1[:, 1, 0]), -(x0[:, 0, 0] - x1[:, 0, 0]),
(x0[:, 0, 0] - x1[:, 0, 0]), (x0[:, 1, 0] - x1[:, 1, 0]),
], axis=1),
[-1, 2, 2],
)
A_I = tf.reshape(
tf.stack([
(x0[:, 1, 0] - x1[:, 1, 0]), (x0[:, 0, 0] - x1[:, 0, 0]),
-(x0[:, 0, 0] - x1[:, 0, 0]), (x0[:, 1, 0] - x1[:, 1, 0]),
], axis=1) / tf.reshape((x0[:, 0, 0] - x1[:, 0, 0])**2 + (x0[:, 1, 0] - x1[:, 1, 0])**2, [-1, 1]),
[-1, 2, 2],
)
c_0 = tf.matmul(A, x0)
c_1 = tf.matmul(A, x1)
c_2 = tf.matmul(A, x2)
c_01 = tf.stack([c_0[:, 1, 0], c_1[:, 1, 0]], axis=1)
c_min = tf.reduce_min(c_01, axis=1)
c_max = tf.reduce_max(c_01, axis=1)
c = tf.reduce_max(
tf.stack([
tf.reduce_min(tf.stack([c_2[:, 1, 0], c_max], axis=1), axis=1),
c_min,
], axis=1),
axis=1,
)
y = tf.reshape(
tf.stack([c_0[:, 0, 0], c], axis=1),
[-1, 2, 1],
)
x4 = tf.matmul(A_I, y)
dist = tf.reduce_sum((tf.reshape(x2, [-1, 2]) - tf.reshape(x4, [-1, 2]))**2, axis=1)
idx = tf.argmin(dist)
res = x4[idx, :, 0]
return res[0].numpy(), res[1].numpy()
# 후보군에 대해 가장 가까운 도로로 재배치 한다.
result = list()
for idx, row in final_gdf.iterrows():
point = [row['geometry'].coords[0][0], row['geometry'].coords[0][1]]
result.append(find_nearest_point(x0, x1, point))
# output list를 df로 변경한다.
result_gdf = gpd.GeoDataFrame(result, columns=['X', 'Y'])
result_gdf['XY'] = result_gdf['X'].apply(lambda x : str(x)) + ' ' + result_gdf['Y'].apply(lambda x : str(x))
result_gdf['geometry'] = result_gdf['XY'].apply(lambda x : Point(float(x.split()[0]), float(x.split()[1])))
result_gdf = result_gdf.set_crs(default_crs)
result_gdf = result_gdf[['geometry']]
result_gdf
# 후보군을 고양시 지도 위에 띄워본다.
from folium import IFrame
from folium.plugins import MarkerCluster
import folium
import random as rd
m = folium.Map(location=center, tiles="OpenStreetMap", zoom_start=12)
width, height = 310,150
popups, locations = [], []
for idx, row in result_gdf.iterrows():
locations.append([row['geometry'].y, row['geometry'].x])
h = folium.FeatureGroup(name='candidate')
h.add_children(MarkerCluster(locations=locations))
m.add_children(h)
# 사용하기 전 변수 세팅
result_gdf = result_gdf.to_crs(meter_crs)
args = [result_gdf, gdf_boundary_town, "DONG_NM"]
# 함수 실행
result_gdf_visual = join_Town2Candidate(*args)
result_gdf_visual.head()
# 행정동별 자전거 정류장 정보를 모두 합해주는 함수
def calculateStationPerTown2(gdf1, gdf2, rel_method=gpd.GeoSeries.intersects):
"""
input : (gdf_station, gdf_pop_dist, latitude, longitude, coverage, space_relation_method, cal_method)
output : ret_gdf
describe : 행정동 데이터 행 하나에 해당하는 정류장의 개수를 feature로 추가한다.
space_relation_method = gpd.GeoSeries.within, gpd.GeoSeries.intersects
"""
ret_df = copy.deepcopy(gdf1)
for i in range(ret_df.shape[0]):
town = gdf1.iloc[i]["geometry"]
temp_df = gdf2[rel_method(gdf2, town)]
if temp_df.shape[0] == 0:
ret_df.loc[i, "STATION_COUNT"] = 0
else:
ret_df.loc[i, "STATION_COUNT"] = temp_df.shape[0]
return ret_df
# 사용하기 전 변수 세팅
args = [gdf_boundary_town, result_gdf_visual, gpd.GeoSeries.intersects]
# 함수 실행
gdf_boundary_town_result = calculateStationPerTown2(*args)
gdf_boundary_town_result.head()
# 후보군을 포함한 행정동 내 자전거 거치대 수를 파악한다.
from branca.colormap import linear
m = folium.Map(center, zoom_start=12,tiles=None)
feature_group2 = folium.FeatureGroup(name='STATION_COUNT',overlay=False).add_to(m)
fs = [feature_group2]#,feature_group1]#,feature_group2]
feature_list = ["STATION_COUNT"]#,"STATION_COUNT"]#,"STATION_PER_POP"]
for i in range(len(feature_list)):
choropleth1 = folium.Choropleth(
geo_data=gdf_boundary_town,
name='choropleth',
data=gdf_boundary_town_result,
columns=['DONG_NM', feature_list[i]],
key_on='feature.properties.DONG_NM',
fill_color='YlGn',
nan_fill_color="black",
fill_opacity=0.7,
line_opacity=0.2,
legend_name='Arrival (in Quintal)',
highlight=True,
line_color='black').geojson.add_to(fs[i])
colormap = linear.YlGn_09.scale(
gdf_boundary_town[feature_list[i]].min(),
gdf_boundary_town[feature_list[i]].max()).to_step(10)
colormap.caption = feature_list[i]
colormap.add_to(m)
folium.TileLayer('cartodbpositron',overlay=True,name="background").add_to(m)
folium.LayerControl('topright',collapsed=False).add_to(m)
m
# 거치대 용량 평균
np.mean(station_df["HOLDER_CAPACITY"])
거치대 용량의 평균이 26인 것을 확인했다.
우리는 이 거치대 용량의 평균에 맞는 숫자를 배정하려 했다.
결과적으로 상대 수요를 기반으로 정렬하고 이를 3등분하여 순차적으로 30, 25, 20개의 거치대 수량을 배정하여 이 평균치에 맞도록 하였다.
# 상대 수요 기반으로 거치대 수량 배정
station_num_df = pd.DataFrame(np.abs(final_data_y), columns=['value'])
station_num_df['station_id'] = station_num_df.index
station_num_df['num'] = pd.qcut(station_num_df['value'], 3, labels=[20, 25, 30])
station_num_df
result_gdf = result_gdf.to_crs(default_crs)
# 답안 제출
result_gdf['스테이션 번호'] = result_gdf.index
result_gdf['거치대 수량'] = station_num_df['num']
result_gdf['X 좌표(위도)'] = result_gdf['geometry'].apply(lambda x : x.coords[0][1])
result_gdf['Y 좌표(경도)'] = result_gdf['geometry'].apply(lambda x : x.coords[0][0])
result_gdf = result_gdf[['스테이션 번호', '거치대 수량', 'X 좌표(위도)', 'Y 좌표(경도)']]
result_gdf
기존에 설치된 거치대를 기반으로 새로운 후보군의 수요를 예측했다.
하지만 기존의 정류장 수 160개를 기반으로 이러한 수요 특성을 파악하기에는 모집단의 수가 너무 적다는 문제가 있다. 최대한 일반화 성능을 높히는 방향으로 진행했지만, 보다 많은 수의 데이터가 있었다면 더 정확한 제안을 하지 않았을까 한다.
시간 복잡도의 문제로, 지역을 기반으로 Rebalancing을 해결하고자 하였다. 우리가 선정한 Metric은 근본적으로 3차원 그래프의 부피를 구하는 작업이다. 이 작업은 굉장히 많은 연산량을 요구한다. 따라서 지역의 Rebalancing을 최소화하는 후보군을 선택하고, 이것들의 모임을 전체 고양시의 거치대의 불균형 문제를 해결하는 답안으로 제출했다. 더 좋은 Metric을 설정하거나 빠른 연산이 가능하다면 전역 Rebalancing을 해결하는 답안을 찾고 싶다.
지역적으로 후보군을 찾는데 있어서 역시 완전 탐색을 통한 최적 후보군을 찾지 못했다. 이 역시 연산량에서 근간한 문제로, 완전 탐색으로 대표 후보군 3N개로부터 N개의 최적 후보군을 찾기 위해서는 $O(n!)$ 정도의 시간 복잡도를 요구한다. 또한 각각의 Metric을 구하는 것도 6.2에서 말했듯 굉장히 높은 연산량을 요구하기 때문에 이를 시도하지 못했다.
이를 보완 할 수 있는 알고리즘으로 유전 알고리즘을 제안한다. 유전 알고리즘은 후보군을 랜덤으로 여러개 생성한뒤 Metric을 통해 세대를 거듭하고 이에 따라 더 낮은 Metric을 갖는 후보군을 찾는다.
최적 입지를 선정하는데 있어서 수요만을 고려한다면, 우리가 진행한 평균 수요를 기반으로 한 방법은 매우 합리적이다. 하지만 Rebalancing 문제에 있어서 평균 수요를 가지고 이를 고려하는 것은 정적인 불균형 문제만을 해소하기 때문에 약간의 부족함이 있다. 만약 시간대 별로 자전거 재배치 차량의 동선, 수량 등에 대한 데이터가 있었다면, 좀더 나은 Rebalancing 해결 방안을 제시할 수 있지 않았을까 하는 아쉬움이 남는다.